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PREFACE 


The  work  reported  herein  represents  the  interim  technical  report  at  the  completion 
of  Phase  I  of  the  Air  Force  Systems  Command,  Aeronautical  Systems  Division  contract 
F33615-86-C-3006  "Low  Density  Real  Gas  Flows  About  Hypersonic  Vehicles."  The  Air 
Force  Project  Engineer  is  Arthur  B.  Lewis.  The  work  was  conducted  by  the 
Computational  Fluid  Dynamics  group  of  the  Boeing  Aerospace  Company  (BAC) 
Engineering  Technology  Organization.  The  BAC  project  manager  is  Dr.  T.  C.  Nark,  and 
the  deputy  program  manager  is  Dr.  S.  F.  Birch.  BAC  technical  leads  for  Phase  1  were 
Dr.  J.  J.  Hoffman  and  R.  G.  Hopcroft  (algorithm)  and  Dr.  T.  R.  Bussing  and  R.  S.  Wong 
(chemistry). 


|  Accesion  For  1 

NTfS 

CRA&I 

a — 

one 

TAB 

□ 

Liidn' 

o.ioced 

□ 

JllStlfi 

:uticri 

By 

Dist  ib 

utio  :  / 

Availability  Codes 

Avd'l 

■rid !  or 

Dist 

Special 

M 

iii/iv 


TABLE  OF  CONTENTS 


Page 


L  INTRODUCTION  .  1 

H.  COMPUTATIONAL  ALGORITHM  .  3 

2.1  INTRODUCTION  .  3 

2.2  NAV1ER-STOKES  SOLUTION  ALGORITHM  .  5 

2.3  NONEQUILIBRIUM  CHEMISTRY  ALGORITHM  .  26 

2.4  VECTORIZATION  AND  PARALLEL  PROCESSING  .  39 

m.  AIR  CHEMISTRY  REACTION  MODEL  AND  WALL  CATALYSIS  .  44 

3.1  INTRODUCTION  .  44 

3.2  DESCRIPTION  OF  GENERALIZED  EQUATIONS  .  45 

3.3  CHEMICAL  REACTION  MODEL  DEVELOPMENT  .  51 

3.4  WALL  CATALYSIS  .  76 

IV.  LEESIDE  EFFECTS  AND  TURBULENCE  MODELS  .  78 

4.1  LEESIDE  EFFECTS  .  78 

4.2  TURBULENCE  MODELS  .  81 

V.  CONCLUSIONS  AND  RECOMMENDATIONS  .  84 

REFERENCES  .  86 

ABBREVIATIONS  .  94 

NOMENCLATURE  .  95 


v 


LIST  OF  ILLUSTRATIONS 

Figure  Page 

1  50  x  20  Grid  in  Two-Dimensional  Plane  for  Three- 

Dimensional  Compression  Corner  Oblique  Shock  Test  Case  21 

2  Calculated  Results  for  Three-Dimensional  Compression 
Corner  Test  Case:  (a)  Mach  Contours,  (b)  Pressure  Contours, 

(c)  Temperature  Contours  22 

3  35  x  30  Grid  in  Two-Dimensional  Plane  for  Three-Dimensional 

Blunt  Nose  Normal  Shock  Test  Case  24 

4a  Mach  Contours  for  Three-Dimensional  Blunt  Nose  Test  Case  24 

4b  Pressure  Contours  for  Three-Dimensional  Blunt  Nose  Test  Case  25 

4c  Temperature  Contours  for  Three-Dimensional  Blunt  Nose 

Test  Case  25 

5  Mass  Concentration  on  Stagnation  Streamline  34 

6a  Two-Dimensional  Contours  of  Mass  Concentration  of  N2  (Percent)  34 

6b  Two-Dimensional  Contours  of  Mass  Concentration  of  O2  (Percent)  35 

6c  Two-Dimensional  Contours  of  Mass  Concentration  of  NO  (Percent)  35 

6d  Two-Dimensional  Contours  of  Mass  Concentration  of  N  (Percent)  36 

6e  Two-Dimensional  Contours  of  Mass  Concentration  of  O  (Percent)  36 

7  Temperature  Ratios  on  Stagnation  Streamline  37 

8  Two-Dimensional  Temperature  Contours  (Kelvin): 

(a)  Translational- Rotational  Temperature,  (b)  Vibrational 
Temperature  of  N2  38 

9  Multiprocessing  the  Explicit  Algorithm  on  a  Cray  X-MP  42 

10  Subtasks  for  Developing  a  Chemistry  Reaction  Model  and 

the  Selection  of  Reaction  Rates  and  Equilibrium  Constants  44 

11  Thermal  and  Chemical  Nonequilibrium  Regions  About 

a  Stagnation  Point  45 

12  Defining  Equations  for  Multi-Temperature  Hypersonic  Flows  50 

13  Equation  Modeling  Hierarchy  50 

14  Summary  of  Air  Chemistry  Model  Development  Process  52 


vi 


LIST  OF  ILLUSTRATIONS  (Continued) 


Figure 

Page 

15 

Comparison  of  KD,  W,  and  PM  Model  Calculations  With 

Peak  Nitric  Oxide  Measurements  of  Camae  and  Feinberg 

60 

16 

Comparison  of  KD,  W,  and  PM  Model  Calculations  With  Camac 
and  Feinberg's  Measured  Nitric  Oxide  "Time-of-Peak"  Periods 

61 

17 

W,  PM,  and  KD  Electron  Density  Calculations  (particles/cm3) 
Versus  Distance  (cm)  Downstream  of  a  Shock  Front  for 

Mach  15.6  Using  Esch  Thermodynamic  Properties 

61 

18 

W,  PM,  and  KD  Electron  Density  Calculations  (particle/cm'*) 
Versus  Distance  (cm)  Downstream  of  a  Shock  Front  for 

Mach  18.5  Using  Esch  Thermodynamic  Properties 

62 

19a 

W,  PM,  and  KD  Electron  Density  Calculations  (particles/cm'*) 
Versus  Distance  (cm)  Downstream  of  a  Shock  Front  for 

Mach  20.2  Using  Esch  Thermodynamic  Properties 

62 

19b 

W  and  PM  Electron  Density  Calculations  (particles/cm3) 

Versus  Distance  (cm)  Downstream  of  a  Shock  Front  for 

Mach  20.2  Using  Esch  Thermodynamic  Properties 

63 

20 

W,  PM,  and  KD  Electron  Density  Calculations  (particles/cm3) 
Versus  Distance  (cm)  Downstream  of  a  Shock  Front  for 

Mach  15.6  Using  Modified  Esch  Thermodynamic  Properties 

65 

21 

W,  PM,  and  KD  Electron  Density  Calculations  (particles  /cm3) 
Versus  Distance  (cm)  Downstream  of  a  Shock  Front  for 

Mach  18.5  Using  Modified  Esch  Thermodynamic  Properties 

65 

22a 

W,  PM,  and  KD  Electron  Density  Calculations  (particles/cm3) 
Versus  Distance  (cm)  Downstream  of  a  Shock  Front  for 

Mach  20.2  Using  Modified  Esch  Thermodynamic  Properties 

66 

22b 

W  and  PM  Electron  Density  Calculations  (particles/cm3) 

Versus  Distance  (cm)  Downstream  of  a  Shock  Front  for 

Mach  20.2  Using  Modified  Esch  Thermodynamic  Properties 

66 

23 

W  Electron  Density  Calculations  (particles/cm3)  Versus 

Distance  (cm)  Downstream  of  a  Shock  Front  for  Mach  15.6, 

18.5,  and  20.2  Using  Modified  Esch  Thermodynamic  Properties 

67 

24 

PM  Electron  Density  Calculations  (particles/cm3)  Versus 

Distance  (cm)  Downstream  of  a  Shock  Front  for  Mach  15.6, 

18.5,  and  20.2  Using  Modified  Esch  Thermodynamic  Properties 

68 

LIST  OP  ILLUSTRATIONS  (Concluded) 


Figure 

Page 

25 

Modified  PM  Electron  Density  Calculations  (particles/cm^) 

Versus  Distance  (cm)  Downstream  of  a  Shock  Front  for 

Mach  15.6,  18.5,  and  20.2  Using  Modified  Esch 

Thermodynamic  Properties 

70 

26 

Typical  Transport  Properties  Methodology 

71 

27 

Specie  Catalysis  at  the  Wall  Surface 

77 

28 

Behavior  of  Leeside  Function  F(y)  at  High  Angle-of-Attack 

79 

29 

Normalized  Pitot  -  Pressure  Contours  in  Flowfield  for  a  Yawed 
5-deg  Cone  at  M®  =  1.8,  «  =  12.5  deg,  Rex  =  28.9  (106)): 

(a)  Experiment,  ref.  52,  (b)  Computed,  ref.  51. 

80 

30 

Computed  Crossflow  Plane  Velocity  Vectors  on  a  5-deg  Cone 
(M.  =  1.8;  «  =  12.5  deg,  Rex  =  28.9(106)) 

80 

viii 


TABLES 


Table  Page 

1  Wray  Air  Model  53 

2  Park  and  Menees  Air  Model  54 

3  Dunn  and  Kang  Air  Model  55 

4  Bittker  Air  Model  56 

5  Bortner  Air  Model  57 

6  Summary  of  Species,  Paths,  and  Rate  Constants 

for  Various  Air  Models  58 

7  Modified  Park  and  Menees  Air  Model  69 

8  Specie  Dynamic  Viscosity  72 

9  Specie  Thermal  Conductivity  72 

10  Binary  Diffusion  Coefficient  of  Specie  O  into  Specie  s  73 

11  Binary  Diffusion  Coefficient  of  Specie  O2  into  Specie  s  73 

12  Binary  Diffusion  Coefficient  of  Specie  NO  into  Specie  s  73 

13  Binary  Diffusion  Coefficient  of  Specie  N  into  Specie  s  73 

14  Binary  Diffusion  Coefficient  of  Specie  NO+  into  Specie  s  73 

15  Binary  Diffusion  Coefficient  of  Specie  N2  into  Specie  s  73 

16  Binary  Diffusion  Coefficient  of  Specie  N+  into  Specie  s  74 

17  Binary  Diffusion  Coefficient  of  Specie  0+  into  Specie  s  74 

18  Candidate  Leeside  Evaluation  Cases  81 


L  INTRODUCTION 


Renewed  interest  in  hypersonic  flow  calculations  after  nearly  20  years  of  reduced 
activity  has  revealed  that  our  current  calculation  capability  is  effectively  limited  to 
inviscid  flows  or  relatively  benign  viscous  flows.  Real  hypersonic  flows,  however,  have 
significant  viscous  and  nonequilibrium  chemistry  effects.  In  many  cases,  the  viscous 
layer  constitutes  a  large  fraction  of  the  shock  layer,  and  chemical  nonequilibrium 
effects  are  strong. 

There  are  currently  no  known  computational  analysis  programs  which  can 
adequately  and  efficiently  describe  the  low  density  real  gas  flows  about  hypersonic 
vehicles.  It  is  generally  accepted  that,  for  continuum  flows,  the  full  Navier-Stokes 
equations  accurately  describe  the  fluid  dynamic  system.  Most  current  computational 
programs,  however,  are  generally  based  on  some  reduced  form  of  the  Navier-Stokes 
equations,  such  as  the  thin  layer  Navier-Stokes  (TLNS),  parabolized  Navier-Stokes 
(PNS),  and  reduced  Navier-Stokes  (RNS)  equations.  Techniques  based  on  these  reduced 
equation  sets  necessarily  contain  simplifying  assumptions  which  may  invalidate  their  use 
for  many  problems  of  current  interest.  Additionally,  many  of  the  techniques  based  on 
the  reduced  equation  sets  are  not  suitable  for  including  reacting  gas  chemistry  effects. 

The  goal  of  the  current  contract  effort  is  directed  at  producing  a  computer  code 
which  can  accurately  and  efficiently  predict  trends  in  heating  rates,  pressures,  forces, 
and  moments  on  complex  shapes  in  low-density  flow  at  high  Mach  numbers,  and  will 
provide  a  realistic  means  for  identifying  real  gas  effects  to  assist  in  optimizing  vehicle 
aerodynamics  for  hypersonic  flight.  This  program  is  based  on  a  coupled,  three- 
dimensional,  full  Navier-Stokes/chemistry  solution  algorithm,  and  will  solve  the  fluid 
dynamic  and  chemistry  fields  about  hypersonic  vehicles  150-200  nose  radii  long,  with  a 
target  computational  time  of  one  hour. 
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The  contract  is  divided  into  two  phases:  (I)  Development  of  a  Flow  Model  and 
Computational  Algorithm,  and  (II)  Program  Development.  The  work  reported  herein 
describes  the  status  of  the  contract  at  the  completion  of  Phase  I.  Phase  I  consisted  of 
six  tasks:  (1)  computational  algorithm,  (2)  leeside  effects,  (3)  development  of  a 
chemical  reaction  model,  (4)  selection  of  reaction  rates  and  equilibrium  constants,  (5) 
turbulence  models,  and  (6)  wall  catalysis.  Phase  I  results  for  Task  (1)  will  be  described 
in  Section  II,  while  results  for  Tasks  (3)  and  (4)  will  be  described  in  Section  III,  and 
results  for  Tasks  (2)  and  (6)  will  be  discussed  in  Section  IV.  Section  V  presents 
conclusions  from  Phase  I  results,  and  recommendations  based  on  these  results. 

The  work  accomplished  in  Phase  I,  and  reported  in  this  technical  report,  forms  a 
strong  basis  for  the  Phase  II  program  development  effort,  and  will  result  in  a  computer 
code  which  will  accurately  and  efficiently  calculate  low  density  real  gas  flows  about 
hypersonic  vehicles. 


H.  COMPUTATIONAL  ALGORITHM 


2.1  INTRODUCTION 

2.1.1  Statement  of  Work 

Contract  Task  4.2.1,  "Computational  Algorithm",  required  development  of  a 
computational  algorithm  for  the  study  of  low  density  real  gas  flows.  The  algorithm 
development  consisted  of  three  components:  (1)  a  solution  algorithm  for  the  fluid 
dynamic  (Navier-Stokes)  equations,  (2)  a  solution  algorithm  for  appropriate  reacting  gas 
chemistry  equations,  and  (3)  vectorization  and  parallel  processing  considerations  for  the 
algorithms  in  (1)  and  (2).  As  stated  by  the  request  for  proposal  (RFP),  the  fluid 
dynamics  of  the  shock  layers  and  leeside  effects  requires  full  diffusion  terms  in  the 
Navier-Stokes  (NS)  equations.  Additionally,  the  real  gas  chemistry  effects  requires 
tracking  of  electron  densities  as  well  as  species  in  the  flowfield.  The  resulting 
complexity  of  the  coupled  NS/chemistry  equations  requires  that  parallel  processing  and 
vectorization  of  the  algorithms  be  utilized  to  improve  the  runtime  of  the  solution 
procedure. 

2.1.2  Goals 

The  goal  of  the  computational  algorithm  task  was  to  develop  NS  and  chemistry 
algorithms  which  could  be  used  to  accurately  and  efficiently  calculate  the  low  density 
real  gas  flowfield  about  hypersonic  vehicles  150  to  200  nose  radii  long.  Accurate 
solution  of  the  fluid  dynamics  of  such  flows  requires  resolution  of  shocks,  boundary 
layers,  and  leeside  effects,  while  accurate  solution  of  chemistry  effects  requires 
resolution  of  species  concentrations  and  electron  densities. 

2.1.3  Approach 

The  approach  selected  for  the  computational  algorithm  task  of  Phase  I  consisted  of 
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developing  the  NS  solution  algorithm  and  the  chemistry  algorithm  separately  (using 
existing  codes),  while  simultaneously  performing  an  analysis  of  parallel  processing  and 
vectorization  requirements  for  both  algorithms.  Prototype  NS  and  chemistry  algorithms 
based  on  the  MacCormack  implicit  algorithm  (Refs.  1,  2)  were  chosen  for  further 
development.  The  selection  of  these  prototype  algorithms  will  be  discussed  in  detail  in 
Sections  2.2  and  2.3. 

Phase  I  development  of  the  NS  algorithm  consisted  of  demonstrating  the  use  of  the 
MacCormack  implicit  algorithm  for  calculating  high  Mach  number  flows.  Chemistry 
algorithm  development  during  Phase  I  consisted  of  evaluating  MacCormack's  (Ref.  2) 
two-dimensional  (2D)  high  mach  number  calculations,  studying  the  extension  of  this 
algorithm  to  three-dimensions,  and  studying  the  coupling  of  the  chemistry  algorithm 
with  the  NS  algorithm.  Vectorization  and  parallel  processing  requirements  were  studied 
throughout  the  course  of  Phase  I,  and  included  evaluations  of  hardware  and  software 
systems.  Hardware  evaluations  of  a  number  of  different  computing  machines  were 
performed,  and  multiprocessing  software  was  also  evaluated.  The  components  of  the 
computational  algorithm  (NS  algorithm,  chemistry  algorithm,  vectorization  and  parallel 
processing  requirements)  are  described  separately  in  Sections  2.2,  2.3,  and  2.4, 
respectively.  The  intent  here  is  to  combine  the  results  of  these  studies  during  the  Phase 
II  program  development  effort. 

2.1.4  Phase  I  Computational  Algorithm  Summary 

The  Phase  I  computational  algorithm  development  resulted  in  a  NS  algorithm 
capable  of  calculating  hypersonic  flows,  and  a  chemistry  algorithm  compatible  with  the 
NS  algorithm  also  capable  of  calculating  hypersonic  flows.  Both  algorithms  have 
calculated  flows  about  a  blunt  body  near  Mach  20.  The  hypersonic  calculations  are 
presented  in  Sections  2.2  and  2.3.  Computer  hardware  and  software  evaluations  in 
Phase  1  assessed  the  state-of-the-art  in  vectorization  and  parallel  processing,  and  parts 
of  the  prototype  NS  algorithm  were  vectorized  and  multiprocessed.  These  results  are 
discussed  in  Section  2.4. 
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2.2  NAVIER-STOKES  SOLUTION  ALGORITHM 


2.2.1  Background 

Accurate  calculation  of  the  reacting  gas  flowfield  about  a  hypervelocity  vehicle  at 
high  altitude  requires  an  accurate  description  of  the  fluid  dynamics  of  the  flowfield. 
The  fluid  dynamics  of  hypersonic  flows  are  characterized  by  thin  shock  layers  which  can 
be  fully  viscous.  The  viscous  region  affects  the  shock  structure,  which  in  turn  affects 
the  inviscid  flow,  leading  to  strong  three-dimensional  (3D)  viscous/inviscid  interactions 
(Refs.  3,  4).  Some  form  of  the  Navier-Stokes  equations  must  be  used  to  model  the  flow 
since  viscous  effects  will  be  important. 

There  are  a  number  of  different  forms  of  the  Navier-Stokes  equations  currently  in 
use.  Four  of  the  most  popular  include  parabolized  Navier-Stokes  (PNS),  reduced  Navier- 
Stokes  (RNS),  thin  layer  Navier-Stokes  (TLNS),  and  full  Navier-Stokes  (FNS).  PNS 
methods  assume  that  upstream  influence  is  negligible,  and  are  therefore  not  suited  for 
the  current  study  since  the  strong  viscous/inviscid  interactions  expected  in  the 
hypersonic  flows  of  interest  will  create  a  strong  streamwise  dependence.  RNS  methods 
attempt  to  add  streamwise  dependence  to  the  basic  PNS  equations  through  a  pressure 
relaxation  technique  (Ref.  5).  These  methods  are  in  general  valid  only  for  flows  with 
weak  streamwise  interaction.  For  the  current  study,  the  extent  of  streamwise 
dependence  may  be  large,  especially  in  separated  regions  on  the  leeside  of  the  body. 
Thus,  RNS  methods  are  also  not  suitable  for  the  current  study.  TLNS  methods  retain 
streamwise  derivatives,  but  only  retain  viscous  terms  in  a  direction  normal  to  a  surface 
(Ref.  6).  For  the  problems  to  be  considered  in  this  contract,  the  three-dimensionality  of 
the  viscous  shock  layer,  as  well  as  the  leeside  separation  region,  requires  that  all  three 
components  of  diffusion  be  modeled.  Thus,  TLNS  methods  will  not  satisfy  the 
requirements  for  this  contract. 

The  primary  motivation  for  using  any  of  the  above  reduced  forms  of  the  Navier- 
Stokes  equations  (PNS,  RNS,  TLNS)  is  the  savings  in  computer  time  and  storage 
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requirements  resulting  from  such  reductions.  These  reductions  are  necessarily 
accompanied  by  restrictions  to  their  applications,  making  the  reduced  forms  of  the 
Navier-Stokes  equations  valid  only  for  a  certain  range  of  problems.  The  current 
contract  objective,  to  calculate  three-dimensional  viscous  flows  with  complex 
viscous/inviscid  interactions,  does  not  completely  satisfy  the  requirements  for  any  of 
the  reduced  forms  described  above.  Thus,  the  full  Navier-Stokes  equations  are  required 
to  satisfy  the  objectives  of  this  contract. 

Until  recently  (in  the  past  10  years),  it  had  been  nearly  impossible  to  numerically 
solve  the  full  Navier-Stokes  equations  for  any  realistic  aerodynamic  configurations,  due 
largely  to  the  prohibitive  computer  time  and  storage  requirements  needed  for  grids 
capable  of  producing  accurate  flowfield  calculations.  With  the  emergence  of 
"supf  .'computers",  it  is  becoming  feasible  to  solve  the  full  Navier-Stokes  equations  for 
realistic  geometries  in  a  cost/time  frame  competitive  with  wind  tunnel  testing. 
Considering  the  fact  that  complex  hypersonic  flows  are  difficult,  if  not  impossible,  to 
model  in  the  wind  tunnel  (and  therefore  very  expensive),  it  can  be  seen  that  numerically 
solving  the  full  Navier-Stokes  equations  for  real  problems  of  engineering  interest  is  not 
as  forbidding  as  it  had  been  a  decade  ago. 

Having  selected  the  appropriate  equation  set  (i.e.  the  full  Navier-Stokes  equations), 
it  was  necessary  to  select  a  solution  technique.  In  the  past  10-15  years,  two  solution 
methods  for  the  full  Navier-Stokes  equations  have  received  most  of  the  computational 
fluid  dynamics  (CFD)  community's  attention:  (1)  the  approximate  factorization 
alternating  direction  implicit  (ADI)  technique  of  Beam  and  Warming  (Refs.  7,  8),  and  (2) 
the  implicit  unfactored  Gauss-Seidel  line  relaxation  technique  of  MacCormack  (Refs. 
1,  9). 

Approximately  factored  (ADI)  techniques  have  long  been  known  to  suffer  a 
factorization  error  which  severely  limits  the  computational  time  step,  while  the 
unfactored  Gauss-Seidel  (GS)  line  relaxation  solution  technique  does  not  possess  such  a 
time  step  limitation.  Additionally,  Yee  and  Shinn  (Refs.  10,  11)  have  noted  that  "there 
appears  to  be  no  straightforward  way  of  utilizing  ADI  approaches  for  nonlinear  system 
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cases  with  general  stiff  source  terms.”  One  of  the  primary  goals  of  the  current  contract 
is  to  calculate  the  reacting  gas  properties  of  the  flowfield,  a  calculation  characterized 
by  the  presence  of  stiff  chemistry  source  terms.  The  MacCormack  unfactored  GS 
algorithm  is  readily  amenable  to  the  incorporation  of  coupled  chemistry  equations  with 
stiff  source  terms,  and  has  already  been  demonstrated  for  2D  flow  (Refs.  2,  10). 

Based  upon  the  above  arguments,  the  MacCormack  implicit  unfactored  GS  line 
relaxation  solution  algorithm  for  the  NS  equations  was  selected  for  development  under 
this  contract. 

2.2.2  Features  of  the  Proposed  Algorithm 

The  MacCormack  implicit  algorithm  for  the  NS  equations  possesses  a  number  of 
features  which  make  it  desirable  for  the  current  calculation  of  iow  density  real  gas 
flows  about  hypersonic  vehicles.  These  features  are  summarized  in  this  section. 

Arbitrarily  Transformed  Coordinate  System.  The  MacCormack  implicit  algorithm 
is  written  in  terms  of  generalized  curvilinear  coordinates,  allowing  body  fitted  meshes 
to  be  used.  Body  fitted  meshes  facilitate  geometry  definition  (even  for  complex  3D 
geometries),  and  also  facilitate  implementation  of  boundary  conditions. 

Arbitrary  Boundary  Condition  Specification.  The  MacCormack  implicit  algorithm 
implements  a  full  range  of  supersonic  boundary  conditions,  including  solid  wall  no-slip, 
solid  wall  free-slip,  symmetry,  inflow,  outflow,  and  freestream  conditions.  These 
boundary  conditions  are  treated  implicitly,  thereby  maintaining  the  stability  and 
robustness  of  the  basic  algorithm. 

Full  Viscous  Terms.  As  mentioned  in  the  RFP,  retention  of  three  components  of 
the  viscous  terms  is  critical  for  accurate  calculation  of  the  viscous  region  in  the  shock 
layer,  as  well  as  for  leeside  effects  where  large-scale  separation  is  expected.  In  the 
MacCormack  implicit  algorithm,  full  viscous  terms  are  retained  in  the  explicit  part  of 
the  algorithm  (including  cross-derivative  terms),  while  the  diagonal  components  of  the 
transformed  diffusion  terms  are  retained  in  the  implicit  part  of  the  algorithm.  The 


7 


combined  algorithm  exhibits  the  character  of  full  viscous  terms,  and  will  allow  accurate 
calculation  of  the  shock  layer  and  leeside  effects. 

Inviscid  Flux  Splitting.  Finite  difference  approximations  to  the  NS  equations  suffer 
from  discretization  errors  which  restrict  the  stability  and  robustness  of  the  formulation 
unless  some  form  of  artificial  dissipation  is  added.  In  the  past,  artificial  dissipation  for 
the  NS  equations  most  usually  consisted  of  some  fourth  order  dissipation  function,  such 
as  that  described  by  Beam  and  Warming  (Ref.  7),  or  a  combination  of  second  and  fourth 
order  dissipation,  such  as  that  proposed  by  Jameson  (Ref.  12).  The  disadvantage  of  such 
artificial  dissipation  functions  is  that  they  are  relatively  empirical,  and  the  constants 
can  be  adjusted  without  regard  to  the  physics  of  the  flowfield. 

Total  Variation  Diminishing  (TVD)  schemes  represent  another  method  for  adding 
artificial  dissipation  to  the  governing  equations.  TVD  schemes  have  been  shown  by  Yee 
(Ref.  13)  to  be  non-oscillatory  near  strong  shocks.  TVD  techniques  are  computationally 
intensive,  though,  and  are  not  feasible  for  use  in  the  present  study  given  the 
computational  time  constraints  specified  by  the  contract. 

More  recently,  there  has  been  renewed  interest  in  the  flux  splitting  procedure  of 
Steger  and  Warming  (Ref.  14).  Initially  developed  for  the  Euler  equations,  flux  splitting 
treats  the  inviscid  terms  as  combinations  of  forward  and  backward  facing  contributions, 
depending  upon  the  sign  of  the  eigenvectors.  In  this  manner,  flux  splitting  closely 
models  the  physics  of  the  flowfield.  Additionally,  since  the  difference  schemes  are  one¬ 
sided,  the  formulation  is  inherently  dissipative  (for  first  order  flux  splitting),  thereby 
eliminating  the  need  for  an  arbitrary  artificial  dissipation  term.  Flux  splitting  the 
inviscid  fluxes  also  increases  diagonal  dominance  of  the  block  tridiagonal  matrix 
obtained  from  the  implicit  governing  equation,  thereby  enhancing  the  stability  of  the 
tridiagonal  matrix  inversion. 

Gauss-Seidel  Line  Relaxation  Solution  Procedure.  By  utilizing  the  unfactored  block 
tridiagonal  matrix  structure  obtained  by  linearizing  the  implicit  governing  equation,  the 
MacCormack  implicit  algorithm  does  not  exhibit  the  time  step  restrictions  found  in  the 
approximate  factorization  techniques.  The  block  tridiagonal  is  solved  using  Gauss- 
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Seidel  (GS)  line  relaxation,  solving  a  tridiagonal  matrix  for  lines  along  a  prescribed 
direction.  GS  sweeps  update  the  solution  in  the  remaining  two  directions.  Convergence 
for  each  time  step  is  obtained  in  approximately  two  to  three  iterations  (sweeps),  due  in 
part  to  the  use  of  inviscid  flux  vector  splitting  which  enhances  the  stability  of  the 
solution  procedure. 

Chemistry  Algorithm  Compatibility.  The  basic  MacCormack  implicit  algorithm  is 
capable  of  incorporating  an  arbitrary  number  of  chemistry  equations.  The  resulting 
fully-coupled  system  more  accurately  represents  the  physics  of  the  system  than  do 
lagged  chemistry  approaches.  Additionally,  the  implicit  treatment  of  the  chemistry 
equations  removes  the  "stiffness"  problem  generally  associated  with  explicit  solutions  to 
the  NS/chemistry  equations.  The  solution  procedure  retains  the  GS  line  relaxation 
technique  described  for  the  NS  equations.  By  generalizing  the  chemistry  equations,  an 
arbitrary  number  of  chemistry  equations  may  be  added  without  affecting  the  basic 
solution  procedure.  Details  of  the  chemistry  algorithm  will  be  discussed  in  Section  2.3. 


2.2.3  Implicit  Algorithm  Development 

The  MacCormack  implicit  algorithm  for  the  Navier-Stokes  equations  without 
chemistry  solves  the  full  set  of  equations  (conservation  of  mass,  x-momentum,  y- 
momentum,  z-momentum,  and  energy)  in  strong  conservation  form: 


al^  +  aP  +  ^G;  +  ^H_,  _o 

dt  8x  dy  dz 


(1) 


where 


IT  =  jp.pu.pv.pw.ej 

represents  the  (conserved)  solution  vector,  with  p  =  density,  u,  v,  w,  =  cartesian 
components  of  velocity,  and  e  =  total  energy  per  unit  volume.  Here,  primes  denote  the 
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nontransformed  (physical)  flux  vectors,  and  superscript  T  denotes  the  transpose  of  the 
solution  vector. 

The  first  step  in  the  development  of  the  implicit  governing  equation  is  to  split  the 
flux  vectors  F,  G',  H*  into  inviscid  and  viscous  contributions: 
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where  the  inviscid  terms  are  contained  in  the  left-hand  vectors,  the  viscous  terms  are 
contained  in  the  right-hand  vectors,  and 


c  —  t 
xy  yx 


i  =  t 
XL  ZX 


X  =  X 

yz  zy 


,(*+£)■ 

\  dy  dx  1 

\  dz  dx  1 

(dv  dw 

dz  dy 


(Id) 


with  p  a  pressure,  T  =  temperature,  p  =  dynamic  viscosity,  pb  =  -<?)p  =  hulk  viscosity, 
and  \  =  coefficient  of  thermal  conductivity.  Pressure  is  related  to  density  and  energy 
through  an  equation  of  state,  p  =  p(p,  e). 

Transformed  Governing  Equations.  For  grids  clustered  near  a  body  the  cartesian 
coordinate  system  in  Equation  (1)  results  in  a  nonuniform  grid,  and  requires  special 
weighting  considerations  in  the  finite  difference  approximations  to  the  governing 
equations.  The  preferred  procedure  is  to  perform  a  generalized  coordinate 
transformation  from  the  cartesian  (x,  y,  z)  coordinate  system  to  a  generalized 
curvilinear  coordinate  system  (£,,  n»  Q.  The  transformed  governing  equations  are  still 
written  in  strong  conservation  form: 


dU  3F  dG  aH 

—  +  —  +  —  +  —  =  0, 

at  a$  an  a< 


with  the  following  definitions: 


(2) 


u  =  irj-1, 


F  = 


r^  +  c^  +  H’^  j~l, 


1 1 


F'nx  +  G'ny  +  H'n. 


(2a) 


G  = 
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H  = 


F\  +  G'^  +  H*^ 


where  $x>  Sy*  £z>  nx.  <ly.  Hz,  <x>  <y,  <^z  are  metrics  of  the  transformation,  J  is  the  Jacobian 
of  the  transformation,  and  the  nontransformed  (primed)  flux  vectors  are  defined  in 
Equation  (1).  The  transformed  coordinate  system  is  uniform,  and  finite  differences  may 
be  applied  directly  to  the  governing  equations  without  weighting  functions. 

Formulation  of  the  Implicit  Governing  Equation.  The  procedure  for  deriving  the 
implicit  governing  equation  from  Equation  (2)  is  summarized  below.  The  derivation  is 
begun  by  taking  the  time  derivative  of  Equation  (2)  and  interchanging  the  order  of 
differentiation.  Next,  Jacobians  of  the  inviscid  flux  vectors  are  taken  with  respect  to 
the  solution  vector  U,  and  the  viscous  flux  vectors  are  written  directly  in  terms  of  aU/at. 
Multiplying  by  At  to  obtain  the  delta  form  of  the  governing  equations  yields: 
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where  A,  B,  C  are  the  Jacobians  of  the  inviscid  flux  vectors  with  respect  to  the  solution 
vector,  F^,  Fn,  F^,  G^,  Gn,  G<,  H^,  Hn,  H<,  are  the  nine  components  of  the  (transformed) 
viscous  matrices,  and  Nf  is  a  matrix  relating  the  nonconservative  solution  vector  V 
(  =  [p,  u,  v,  w,  ei]T)  to  the  conservative  solution  vector  U.  Note  that  the  internal  energy 
(e i )  per  unit  mass  is  related  to  the  total  energy  per  unit  volume  (e) 
by  e  =  p[ej  +  i(u2  +  y2  +  w2)]. 
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Expanding  the  time  derivative  in  Equation  (3)  in  a  forward  difference  yields: 
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Defining  the  following: 


(4) 


AU"  -  explicit  change  = 


and  5Un  +  l  =  implicit  change  = 


and  substituting  Equations  (4)  -  (5)  into  Equation  (3)  yields  (after  rearranging): 
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where  I  is  the  identity  matrix.  Equation  (6)  represents  the  fully-implicit  governing 
equation  for  the  3D  NS  equations. 

Flux  Splitting  Add  -on.  As  was  discussed  earlier,  most  finite  difference 
representations  of  the  NS  equations  require  some  form  of  artificial  dissipation  to 
maintain  numerical  stability.  Inviscid  flux  vector  splitting  has  been  selected  to  provide 
the  artificial  dissipation,  since  splitting  more  closely  models  the  flow  physics  than 
empirical  artificial  dissipation  functions. 

Flux  vector  splitting  is  accomplished  at  the  computational  cell  face  level.  The  flux 
at  each  cell  face  is  comprised  of  two  components,  one  forward  facing  (associated  with 
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positive  eigenvalues),  and  one  backward  facing  (associated  with  negative  eigenvalues). 
Non-flux-split  procedures  approximate  the  flux  at  a  cell  face  as  a  single  component, 
effectively  central  differencing  the  flux.  First  order  flux  splitting  schemes  forward 
difference  negative  eigenvalue  flux  contributions  and  backward  difference  positive 
eigenvalue  flux  contributions.  For  subsonic  flows  where  both  positive  and  negative 
eigenvectors  exist,  the  procedure  is  equivalent  to  the  central  difference  form.  For 
supersonic  flows,  where  the  eigenvectors  are  all  the  same  sign,  the  flux  splitting  results 
in  a  one-sided  difference.  This  one-sided  difference  will  more  effectively  capture  the 
physics  of  shock  waves  (where  discontinuities  in  properties  exist),  while  at  the  same 
time  providing  the  necessary  dissipation  to  stabilize  the  solution  procedure. 

Flux  splitting  is  implemented  into  Equation  (6)  by  first  diagonalizing  the  inviscid 
flux  vectors  (A,  B,  C)  and  splitting  the  eigenvectors  into  positive  and  negative 
contributions: 


a=saa,si'  =  sa(a;  +a;)sa_,  =  a"  +  a_’ 
b  =  sbabV=Sb(ab*  ^  ab  )  Sb‘  =  B*  +  B- .  (7) 

C  =  SC  ACSC  1  =  Sc(AC  +  AC  )  Sc’  =  C*  +  C‘  • 

where  S,  S~1  represent  the  matrices  which  diagonalize  the  inviscid  flux  vectors,  A 
represents  the  diagonalized  matrix  of  eigenvalues,  and  superscripts  +  and  -  indicate 
positive  and  negative  contributions,  respectively.  The  significance  of  this  splitting 
operation  will  become  apparent  when  the  governing  equations  are  finite  differenced,  as 
the  positive  contributions  will  be  backward  differenced,  and  the  negative  contributions 
will  be  forward  differenced. 

Substituting  the  flux  split  inviscid  terms  in  Equation  (7)  into  the  implicit  governing 
Equation  (6)  yields  the  flux  split  form  of  the  governing  equations: 
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where  the  first  set  of  brackets  surround  the  inviscid  terms,  and  the  second  set  of 
brackets  surround  the  viscous  terms.  Here,  8Un+l  represents  the  implicit  change,  and 
AUn  represents  the  explicit  change,  as  defined  in  Equation  (5). 

Solution  Procedure.  The  first  step  in  the  solution  process  consists  of  discretizing 
the  governing  equation  (Equation  (8))  by  approximating  the  derivative  terms  with  finite 
differences.  As  noted  earlier,  the  following  rules  apply  to  the  discretization  process: 
(1)  positive  inviscid  fluxes  are  backward  differenced,  (2)  negative  inviscid  fluxes  are 
forward  differenced,  and  (3)  viscous  fluxes  are  central  differenced.  Using  these  rules, 
the  resulting  finite  difference  approximation  to  Equation  (8)  becomes: 


1 1  +  At  D  A+  +  D  A~  +  D  B+  +  D  B"  +  D  ,C+  +D  ,C‘  5Un+l 
[  -t  +  t  -n  +  n  -<  +< 

+  4t|D-l(FlD.<  +  F,D.,  +  F<D*<H 
+  D-,(GsD.l  +  G,D.„  +  G<D.<b 

+  D-t(»tD.S  +  H„E.,  +  HiD«)Nf|8lJ’*'  =  (9> 
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where  D-  indicates  backward  differences,  D+  indicates  forward  differences,  and 
subscripts  q,  <  indicate  the  difference  direction.  Note  that  D+  D-  is  equivalent  to  a 
central  difference. 

Applying  the  difference  operators  to  Equation  (9)  and  collecting  like  terms  yields 
the  discretized  form  of  the  governing  equations: 


.ii  +  l  t,,n  +  l  ,  „n+-lKITn-+l  .n  +  l  t„n+l 

a,  oU  .  +•  a.  oU.  .  .  +  a.  6U  ,  . 

1  i.j  +  l,k  0  i,j,k  2  i . j  —  1,  k 


(10) 


+  a"+1SU  ,  ,+a^‘SU  ,  ,  +  a"  +  18U .  .  .  _ ,+  a”  +  18U  ,  ,=  AUn  ,  , 

3  i  +  l.J.k  4  i-l.J.k  5  i.J.k  +  1  6  i.J.k-t  i.J.k 

where  the  constants  ai  through  a$  are  defined  as  follows: 


ao  + '  =  l  +  ^ 


A,  .+ 

i-i.l,  k 

(r‘N'U, 

!FlNr) 

+  t.j.  k 

+ 

+ 

B.V,k  + 

(°,Nr 

k 

c, * 

+ 

c.:J,k+i  + 

Kn< 

l.j,k  +  i 

a"  + 1  =  At 


1“  _  .  +  (g  N.  ) 
+  k  V  n  f/ 


i.J  +  t,k 


n  +  l 


a"  +  1=  At 


-  B +  .  .  +  fc  N  ) 

M-t.k  \  n  r/,  , 


n  +  1 


a3  +  I  =  At 


A.  ...+ 

i  +  iJ.k 


(*M 


i  +  *.J.k 


n  + 1 


a"  +  1  -  At 

4 


—  A+ 


*.i.k 


(F!N. 


i-i.j,  k 


n  +•  1 
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with  subscript  1/2  referring  to  cell  faces.  Note  that  only  diagonal  components  of  the 
viscous  fluxes  (F^,  Gn,  H^)  have  been  retained  in  Equation  (10).  There  is  no  restriction  in 
principle  requiring  exclusion  of  the  cross-derivative  terms;  however,  the  implementation 
of  these  terms  is  difficult  due  to  the  coordinate  transformation  metrics.  It  is  not 
expected  that  neglecting  the  cross-derivative  terms  on  the  implicit  side  of  Equation  (10) 
will  significantly  affect  the  stability  of  the  implicit  algorithm,  since  full  viscous  terms 
are  still  retained  in  the  calculation  of  the  explicit  change  (AUn)  in  Equation  (10). 

Equation  (10)  is  nonlinear  due  to  the  dependence  of  the  coefficients  ao  -  ag  on  the 
solution  vector.  Solution  of  Equation  (10)  is  not  practical  by  any  current  solution 
technique.  Two  assumptions  are  needed  to  make  the  solution  feasible.  First,  the 
coefficients  ao  through  ae  in  Equation  (10)  are  linearized  to  time  level  n.  The  resulting 
equation  represents  a  block  septadiagonal,  whose  solution  is  still  not  practical  for  the 
problems  to  be  considered  in  this  contract  using  current  solution  techniques.  Second, 
selected  terms  are  lagged  such  that  the  resulting  equation  reduces  to  a  block  tridiagonal 
matrix,  for  which  efficient  solution  techniques  do  exist. 

There  are  three  choices  for  the  tridiagonal  structure,  corresponding  to  the  three 
coordinate  directions  (denoted  here  as  i,  j,  k).  It  is  desired  to  utilize  the  tridiagonal 
(implicit)  solution  update  along  lines  normal  to  the  steepest  flow  gradients,  such  as  will 
be  found  normal  to  a  solid  wall  and  across  a  shock  wave.  For  simplicity,  it  is  assumed 
here  that  the  direction  of  steepest  gradients  is  normal  to  the  j  =  constant  lines;  the 
generalization  to  the  other  two  directions  (i,  k)  is  straight  forward. 

To  create  a  tridiagonal  structure  in  the  j-direction,  only  terms  with  subscripts 
(i,  j+1,  k),  (i,  j,  k),  and  (i,  j-1,  k)  are  retained  on  the  left-hand  side  of  Equation  (10). 
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Noting  also  that  the  coefficients  ao  through  ag  have  been  linearized  to  time  level  n, 
Equation  (10)  becomes: 


a"8Un+j,  .  +  a"  8Un  +  ’  +  a"8U' 

1  i.j-t-l,  k  0  i,j,k  2  i,j 


n  + 1  _ 


- 1,  k 


AUn  .  - 

‘.J.  k 


a36ur.,i.k 


a"8Un  .  . . 

4  i  -  l.j,  k 


a"  8Un  .  .  .  —  a«  8Un  .  , 

5  i,j,  k  +  I  6  l. j.  k  —  1 


(ID 


Similar  equations  may  be  written  in  the  i  and  k  directions. 

Equation  (11)  represents  a  block  tridiagonal  algorithm  which  is  implicit  in  the  j- 
direction  for  updating  the  solution  vector  U.  Due  to  lagging  the  i-  and  k-direction 
terms,  a  Gauss-Seidel  line  relaxation  solution  procedure  is  used,  with  each  time  step 
consisting  of  two  or  more  sweeps  in  the  i-  and  k-directions,  with  a  tridiagonal  matrix 
solution  for  each  j-line.  Experience  has  shown  that  two  to  three  sweeps  are  sufficient 
to  obtain  a  converged  solution  for  each  time  step. 

MacCormack's  baseline  implicit  algorithm  has  a  predictor-corrector  format  (Ref. 
1),  with  the  following  six-step  solution  process: 
a.  Calculate  the  explicit  delta  (predictor): 


AUn  =  -At 

i.J,  k 


D  .  F  + 


f  D  F  +  D  ,  F  . 

inv  -*-£  tnv  ±t,  vis 


-At 


D  G+  +  D  G~  +  D,  G 

-q  inv  +q  inv  ±r\  vis 


-At  Id 
i 


.  H  +  +  D  H”  +  D,  .  H 

mv  -**4  inv  ±c  vis 


where  D-  indicates  backward  differences,  D+  indicates  forward  differences,  and  D± 
indicates  either  a  forward  or  backward  difference.  Here,  subscript  inv  denotes  inviscid 
terms  and  subscript  vis  denotes  viscous  terms.  The  differencing  scheme  for  the  viscous 
terms  is  alternated  from  the  predictor  to  the  corrector  step,  resulting  in  (effectively) 
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central  differenced  viscous  terms.  Full  viscous  terms  are  retained  in  the  viscous  flux 
vectors.  Note  that  flux  splitting  has  been  applied  to  the  inviscid  flux  terms  in  the 
calculation  of  the  explicit  delta. 

b.  Calculate  the  implicit  delta  (predictor),  by  solving  Equation  (11)  for  6U  n^,  where 
the  overbar  denotes  an  intermediate  solution. 


c.  Intermediate  (predictor)  update: 

8Un  +  f  =  U"  +  8Un. 

i,j,  k  ijjc  lj.k 

d.  Calculate  explicit  delta  (corrector),  by  repeating  step  (a)  with  alternated  difference 
directions. 

e.  Calculate  the  implicit  delta  (corrector),  by  solving  Equation  (11)  for  5U-n^,)  . 

f.  Solution  update: 


un+i  =  * 

i,j,k  a 


Un.,+8Un'fl  +  8UnV 

‘J-*  i,j,k  1J,k 


As  part  of  this  contract  effort,  a  one-step  version  of  the  solution  algorithm  is  being 
investigated.  The  primary  advantage  of  a  one-step  solution  process  is  the  50%  savings 
in  the  computational  time  associated  with  the  predictor/corrector  solution.  The  one- 
step  solution  procedure  is  summarized  below: 
a.  Calculate  explicit  delta: 


AU".  =  -  At 

i.J.k 


D  F+  +  D  .  F“  +  D , .  F 

-t  inv  +  {,  inv  vis 


-At 


D  G.+  +  D,  G"  +  D .  G 

-q  inv  +  i)  inv  ±i)  vis 


-At 


D  .  H.+  +  Dx/H"  +  D^H 

inv  +  {  inv  ±C  vi  8 
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b. 


Calculate  implicit  delta  for  5U 


n+1 

i,j,k 


from  Equation  (11). 


c.  Update  the  solution: 


Un  +  1  =  Un 

i.J.k 


U.k 


+  8U 


n  +  l 


Although  the  one-step  procedure  described  above  appears  to  be  a  logical  choice,  a 
number  of  issues  regarding  the  correct  form  of  the  differencing  schemes  used  in  the 
one-step  procedure  will  need  to  be  addressed.  In  particular,  it  is  not  clear  if  the  central 
differencing  of  the  viscous  terms,  which  was  achieved  by  alternated  one-sided 
differences  in  the  predictor  and  corrector  steps,  must  be  accomplished  in  one  step,  or  if 
the  central  difference  which  occurs  after  every  two  computational  steps  (retaining  the 
predictor-corrector  format)  is  sufficient.  Additionally,  the  similarity  matrices  used  to 
diagonalize  the  inviscid  fluxes  must  be  evaluated  at  face  centers.  In  the  predictor- 
corrector  solution,  these  centered  values  are  obtained  by  alternated  one-sided 
approximations  in  the  predictor  and  corrector  steps.  Again,  the  one-step  procedure  will 
only  obtain  a  centered  value  after  two  computational  steps  if  the  original  predictor- 
corrector  format  is  retained.  These  topics  will  be  addressed  during  Phase  II  of  this 
contract. 


2.2.6  Test  Cases 

The  flow  physics  associated  with  hypersonic  low  density  real  gas  flows  can  be 
broadly  categorized  into  four  parts:  (1)  weak  (oblique)  shocks,  (2)  strong  (normal) 
shocks,  (3)  boundary  layer/heat  transfer,  and  (4)  chemical  nonequilibrium  effects.  Items 
(l)-(3)  pertain  to  the  fluid  dynamics  of  the  flow,  and  are  critical  features  in  the 
determination  of  item  (4).  Two  test  cases  were  created  in  Phase  I  to  evaluate  the 
prototype  Navier-Stokes  (fluid  dynamics)  solution  algorithm.  These  cases  were  selected 
to  attempt  to  isolate  flow  features  (1)  and  (2)  described  above.  Test  Case  1  consisted  of 
a  compression  corner,  and  attempted  to  isolate  a  weak  (oblique)  shock,  while  Test  Case 
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2  consisted  of  a  3D/axisym metric  blunt  nose  region  which  attempted  to  isolate  a  strong 
(normal)  shock.  Emphasis  at  this  stage  of  the  contract  was  on  determining  whether  the 
prototype  NS  solution  algorithm  could  accurately  predict  gross  flow  features  at  high 
Mach  numbers.  No  attempt  was  made  to  analyze  fine  flow  details,  and  for  this  reason 
no  boundary  layer/heat  transfer  studies  were  conducted  at  this  time;  nor  were  there  any 
attempts  to  compare  the  results  obtained  from  the  two  test  cases  to  experimental  data. 
It  was  considered  sufficient  for  present  purposes  to  demonstrate  global  features  of  the 
solver  for  high  Mach  number  flows. 

Teat  Case  It  Compression  Comer.  The  geometry  for  the  3D  compression  corner 
consisted  of  a  1  inch  (0.0254m)  long  flat  plate  followed  by  an  11.3  degree  half-angle 
ramp.  A  coarse  50  x  20  mesh,  shown  in  Figure  1,  was  used  in  the  two-dimensional  plane, 
and  3  grid  cells  were  used  in  the  spanwise  direction.  The  plate  width  was  taken  as  0.1 
inch.  Test  conditions  for  this  test  case  are  summarized  below. 


Mach  number: 

20.0 

Freestream  velocity: 

5906  m/s 

(19377  ft/s) 

Static  pressure: 

40000  N/m2 

(8351bf/ft2) 

Static  temperature: 

216.7K 

(390°R) 

The  Reynolds  number  for  this  case  was  6.2  (106),  based  on  the  plate  length. 

The  test  case  assumed  an  ideal  gas,  and  the  wall  was  considered  to  be  adiabatic. 

v  i 


Figure  1.  50  x  20  Grid  in  Two-Dimensional  Plane  for  Three-Dimensional  Compression 
Corner  Oblique  Shock  Test  Case 


The  solution  converged  in  approximately  80  implicit  steps.  Parts  (a),  (b),  and  (c)  of 
Figure  2  show  the  Mach,  pressure,  and  temperature  contours,  respectively.  The  y-scale 
has  been  expanded  in  Figure  2  for  clarity.  Note  that  the  calculation  predicts  the  gross 


X/L 


(b)  Pressure  Contours 


X/L 


(c)  Temperature  Contours 

Figure  2.  Calculated  Results  for  Three-Dimensional  Compression  Corner  Test  Case 


physical  features  of  the  flowfield,  including  the  leading  edge  shock  which  impinges  the 
shock  from  the  compression  corner.  The  shock  angle  of  15°  compares  well  with  that  for 
an  11.3  degree  half-angle  wedge  given  by  Ref.  IS. 


Test  Case  2:  3D  Blunt  Nose.  The  second  test  case  was  extracted  from  a  geometry 
created  for  the  1987  Parabolized  Navier-Stokes  (PNS)  worksnop,  held  23-24  September 
1987  at  Wright  Patterson  AFB.  The  geometry  used  here  considered  only  the  nose  region 
of  the  PNS  workshop  model,  including  the  1  inch  (0.0254m)  diameter  blunt  nose  and  a  3 
nose  radii  long  section  of  the  5  degree  ramp  attached  to  the  nose.  Again,  the  intent 
here  was  to  study  gross  physical  characteristics  of  an  isolated  strong  (normal)  shock  to 
determine  the  ability  of  the  prototype  algorithm  to  calculate  strong  shocks  at  high  Mach 
numbers.  The  test  conditions  were  the  same  as  for  the  PNS  workshop,  with  the 
exception  that  the  solid  wall  was  assumed  to  be  adiabatic  for  this  case.  The  test 
conditions  are  summarized  below; 

Mach  number: 

Freestream  velocity; 

Static  pressure: 

Static  temperature: 

The  Reynolds  number  (based  on  nose  radius)  was  10*. 

A  35x30  mesh,  shown  in  Figure  3,  was  used  in  the  2D  plane.  Three  cells  covering  15 
degrees  were  used  in  the  circumferential  direction.  The  outer  boundary  was  located 
0.25  nose  radii  upstream  from  the  nose  tip,  and  the  shock  was  captured  by  the 
calculation.  A  converged  solution  was  obtained  in  approximately  100  implicit  steps,  and 
the  results  are  shown  in  Figure  4.  Part  a  of  Figure  4  shows  the  Mach  contours  for  the 
flowfield,  and  clearly  shows  a  well-defined  normal  shock  near  the  nose.  The  shock 
standoff  distance  was  approximately  0.20  times  the  nose  radius,  consistent  with  that 
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6596  m/s  (21633  ft/s) 

79.779N/m2  (1.671bf/ft2) 

270.65K  (488°R) 
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given  by  Liepmann  and  Roshko  (Ref.  16).  Parts  (b)  and  (c)  of  Figure  4  show  pressure  and 
temperature  contours,  respectively. 

2.3  NONEQUILIBRIUM  CHEMISTRY  ALGORITHM 

2.3.1  Background 

Hypersonic  vehicles  at  high  altitude  encounter  an  environment  characterized  by  low 
density  effects  and  reacting  gas  chemistry.  Therefore,  it  is  required  that  a  reacting  gas 
chemistry  equation  set  be  solved.  There  are  four  basic  methods  for  including  the 
chemistry  equations  with  the  fluid  dynamic  (NS)  equations:  (1)  uncoupled  explicit,  (2) 
uncoupled  implicit,  (3)  coupled  point-implicit,  (4)  coupled  fully-implicit.  Of  these  four 
choices  (1)  is  the  simplest  to  implement  but  represents  the  least  amount  of  flow  physics, 
while  (4)  is  the  most  difficult  to  implement,  and  best  describes  the  flow  physics 
associated  with  the  reacting  gas  flowfield.  For  relatively  benign  environments,  the 
lagged  methods  (1),  (2)  are  probably  sufficient;  however,  for  the  hostile  environments  to 
be  encountered  by  high  altitude  hypersonic  vehicles,  it  is  believed  that  it  will  be 
necessary  to  couple  the  chemistry  equation  set  to  the  fluid  dynamic  equation  set. 

Fully-implicit  coupling  of  the  NS/chemistry  equations  has  been  chosen  for 
development  under  this  contract  since,  at  this  time,  it  appears  to  represent  the  best 
prospect  for  alleviating  the  stiffness  problem.  A  system  of  equations  is  described  as 
stiff  if  some  components  of  the  solution  respond  very  quickly  to  system  perturbations 
(i.e.,  have  short  time  scales)  while  others  respond  slowly  (i.e.,  have  long  time  scales). 
The  differences  in  response  times  can  be  due  to  differences  in  response  times  between 
the  various  equation  sets  describing  the  flow,  such  as  the  fluid  dynamic  and  chemistry 
equations,  or  it  can  be  due  to  large  differences  in  the  time  scales  characterizing 
different  regions  of  the  flowfield.  Explicit  solution  techniques  generally  have  stability 
restrictions  that  relate  the  maximum  allowable  time  step  to  the  fastest  characteristic 
time  scale,  and  therefore  often  require  large  numbers  of  computational  steps,  with 
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correspondingly  large  computational  costs.  Implicit  solution  techniques  do  not  exhibit 
these  stability  restrictions  and  can  therefore  use  much  larger  time  steps. 

One  disadvantage  of  using  a  fully-coupled  implicit  approach  is  associated  with  the 
size  of  the  block  matrix  to  be  solved.  For  a  large  number  of  chemical  equations,  the 
additional  cost  to  solve  the  block  matrix  may  be  prohibitive.  For  the  equation  set  to  be 
used  here,  it  is  not  felt  that  the  additional  computational  time  will  be  excessive. 
Development  of  the  coupled  equation  set  is  described  in  Section  III. 

For  the  test  cases  to  be  considered  in  the  present  contract,  nine  species  have  been 
identified  as  critical  for  accurate  flowfield  calculation,  including  electron  density.  The 
fully-coupled  3D  NS/chemistry  equation  set  represents  a  block  14*14  matrix  (original 
5x5  for  NS  plus  9x9  for  chemistry).  Since  the  equation  set  is  fully-coupled,  it  is  desired 
that  the  solution  technique  for  the  coupled  set  be  the  same  as  for  the  fluid  dynamic 
equations.  Thus,  the  MacCormack  algorithm  (Ref.  2),  extended  to  include  an  arbitrary 
number  of  chemical  equations,  has  been  selected  as  the  algorithm  for  solving  the 
coupled  NS/chemistry  equation  set. 

It  should  be  noted  here  that  not  all  nine  species  will  have  significant  effects  in  all 
areas  of  the  flowfield.  Therefore,  an  equation-adaptive  capability  is  being  investigated 
which  would  allow  a  reduced  equation  set  to  be  solved  in  flowfield  regions  where  some 
species  may  not  significantly  affect  the  solution.  In  such  cases,  the  flexibility  of  the 
MacCormack  coupled  NS-chemistry  algorithm  could  be  exploited  by  reducing  the  size  of 
the  block  matrix  to  be  solved,  thereby  resulting  in  a  savings  in  computer  time.  Inclusion 
of  this  feature  will  be  contingent  upon  determination  of  the  generality  of  such  a 
feature,  as  well  as  the  determination  of  whether  the  code  can  automatically  adjust  the 
equation  set,  or  if  the  user  will  be  required  to  specify  the  range  of  validity  for  each 
chemical  species.  It  is  felt  at  present  that  the  equation  set  should  not  be  arbitrarily 
user-specified,  and  that  the  full  equation  set  should  be  retained  unless  some  well- 
defined  general  rules,  such  as  concentration  limits  for  each  species,  can  be  defined  for 
automatically  adjusting  the  equation  set. 


2.3.2  Features  of  the  Solution  Algorithm 

There  are  a  number  of  features  associated  with  extending  the  MacCormack 
algorithm  to  include  chemistry.  Implementation  of  the  fully-coupled  chemistry 
algorithm  is  straightforward  since  it  follows  the  same  procedure  used  for  developing  the 
MacCormack  implicit  algorithm  for  the  NS  equations.  Additionally,  the  fully-coupled 
chemistry  along  each  Gauss-Seidel  solution  line  provides  the  most  accurate  description 
of  the  flow  physics  by  solving  the  fluid  dynamic  and  chemistry  effects  simultaneously. 
As  noted  earlier,  a  third  advantage  (not  unique  to  this  approach)  is  that  the  implicit 
treatment  of  the  chemistry  source  terms  removes  the  chemistry  stiffness  problem,  and 
permits  larger  time  steps  more  consistent  with  an  implicit  solution  algorithm  to  be 
taken.  The  MacCormack  algorithm  extension  to  include  chemistry  is  general,  implying 
that  an  arbitrary  number  of  species  chemistry  equations  may  be  added,  with  the  only 
limitation  being  computational  (i.e.  solving  the  larger  matrix)  rather  than  physical  (i.e. 
any  number  of  species  may  be  modeled).  Finally,  MacCormack  and  Candler  (Ref.  2) 
have  demonstrated  the  chemistry  algorithm  for  a  two-dimensional  blunt  nose  at 
hypersonic  velocity.  Recently,  Shinn,  Yee,  and  Uenishi  demonstrated  this  algorithm  in 
3D  (Ref.  10)  using  a  TVD  approach.  The  MacCormack-Candler  test  case  will  be 
described  in  Section  2.3.4. 

2.3.3  Chemistry  Algorithm  Development 

The  chemistry  algorithm  will  be  discussed  here  in  terms  of  the  coupled  fluid 
dynamics/chemistry  algorithm  for  completeness.  The  MacCormack  algorithm  for  the 
NS  equations  with  coupled  nonequilibrium  chemistry  solves  the  following  equation  set  (in 
strong  conservation  form): 


<3U  8F  8G  3H 

c  c  c  c  • 

—  +  —  +  —  +  — £  =w, 

3t  3x  3y  3z 


(12) 


where  primes  again  denote  nontransformed  (cartesian)  quantities,  and  subscript  c 
distinguishes  the  coupled  fluid  dynamic/chemistry  vectors  in  Equation  (12)  from  the 
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fluid  dynamic  vectors  in  Equation  (1).  Here,  U'e  =  [p,,  ...»  pns,  pu,  pv,  pw,  e]T  denotes  the 

(conservative)  solution  vector,  with  p,,  ...,  pns  representing  the  species  densities  for 

Ns 

species  1  through  Ns  (Ns  =  number  of  species),  p  =  represents  the  total  density, 

j=i 

u,  v,  w  represent  the  Cartesian  velocity  components,  and  e  represents  the  energy  per 
unit  volume.  The  chemistry  source  term  (  W)  is  written  in  the  transposed  form 
W  =  [cbj,  ...»  wNs>  0,  0,  0,  0]T,  where  Wj  -*  d>Ng  represent  the  source  terms  for  each 
species.  The  determination  of  these  source  terms  is  deferred  until  Section  III. 

The  coupled  fluid  dynamics/chemistry  flux  vectors  (F,e,  G'c>  Ifc)  are  again  split 
into  inviscid  and  viscous  parts,  and  are  written  in  cartesian  coordinates  as  follows: 
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where  inviscid  terms  are  contained  in  the  left-hand  vectors,  viscous  terms  are  contained 
in  the  right-hand  vectors,  and 
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with  p  =  pressure,  T  =  temperature,  A  =  coefficient  of  thermal  conductivity,  p  =  dynamic 
viscosity,  pB  =  -ip  =  bulk  viscosity,  hs  =  specific  enthalpy  of  species  s,  and  flj-^UNg, 
v2-»vns,  wj-*wns  represent  the  Cartesian  diffusion  velocities  for  species  1  through  Ns. 
From  Candler  and  MacCormack  (Ref.  2),  it  is  noted  that  the  diffusion  velocities  are 
proportional  to  the  gradients  of  the  species  concentrations,  assuming  that  pressure, 
thermal,  and  body  force  effects  on  the  diffusion  velocities  are  negligible.  Thus,  the 
diffusion  velocities  may  be  written  in  the  form: 
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where  Ds  is  the  specie  concentration  diffusion  coefficient,  and  Cs  is  the  specie 
concentration.  The  translational,  vibrational,  and  rotational  energy  modes  have  been 
omitted  in  the  coupled  governing  equations  described  above.  These  three  energy  modes 
are  assumed  to  be  in  equilibrium,  for  reasons  to  be  discussed  in  Section  III. 

Note  that  a  species  continuity  equation  is  written  for  each  species,  with  no  global 
species  continuity  equation,  as  opposed  to  the  alternative  form  of  writing  Ns~l  species 
continuity  equations  plus  a  global  continuity  equation.  Yee  and  Shinn  (Ref.  11)  have 
noted  that  the  present  form  is  equivalent  to  the  form  utilizing  the  global  species 
continuity  equation,  but  has  the  advantage  that  the  eigenvectors  of  this  system  are 
easier  to  determine. 

Comparing  Equations  (1)  and  (12),  it  can  be  seen  that  the  chemistry  equations  modify 
both  the  inviscid  and  viscous  flux  vectors.  The  global  continuity  equation  in  the  inviscid 
matrices  of  Equation  (1)  is  replaced  with  specie  continuity  equations  for  each  of  the 
species.  In  the  viscous  matrices,  the  continuity  contribution  (which  was  originally  zero)  is 
replaced  by  the  diffusion  velocity  contributions  from  each  of  the  Ng  species.  The  diffusion 
velocities  also  contribute  to  the  energy  term  in  the  viscous  matrices.  Since  the  original  form  of 
the  flux  vectors  is  not  changed,  the  Gauss-Seidel  relaxation  solution  procedure  developed  for 
the  NS  equations  will  also  apply  to  the  coupled  equation  set  in  Equation  (12). 

Solution  Procedure.  The  mathematical  treatment  of  the  coupled  fluid 
dynamics/chemistry  equation  set  is  similar  to  that  described  in  Section  2.2,  and  will  only 
be  summarized  here.  The  cartesian  form  of  Equation  (12)  is  transformed  to  generalized 
curvilinear  coordinates  (1,  n>  c),  and  the  implicit  governing  equation  is  formulated  from 
the  transformed  equation  set.  Once  again,  the  inviscid  flux  vectors  are  split  into 
positive  and  negative  contributions  (based  on  eigenvalue  sign),  so  that  forward  finite 
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differences  can  be  applied  to  negative  contributions  and  backward  finite  differences  can 
be  applied  to  positive  contributions.  The  implicit  governing  equation  for  the  3D  flux- 
split  fully-coupled  NS/chemistry  equation  set  becomes: 
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where  At  is  the  time  step,  1  is  the  identity  matrix,  D-  indicates  backward  differences, 
D+  indicates  forward  differences,  and  Nc  denotes  the  matrix  relating  the  nonconserved 
solution  vector  to  the  conserved  solution  vector.  Note  that,  as  was  the  case  in  the  fluid 
dynamic  algorithm,  only  diagonal  components  of  the  viscous  stresses  are  retained  in  the 
implicit  formulation,  while  full  viscous  terms  are  retained  in  the  explicit  change  (AUen), 
and  subscript  c  denotes  vectors  associated  with  coupled  NS/chemistry  equations. 
Equation  (13)  is  solved  by  the  same  procedure  described  in  Section  2.2.5,  the  only 
difference  being  the  size  of  the  block  tridiagona'  to  be  solved.  The  solution  of  Equation 
(13)  will  provide  accurate  flowfield  calculations  for  low  density  real  gas  flows  about 
hypersonic  vehicles. 


2.3.4  Test  Case 

The  coupled  NS/chemistry  algorithm  has  been  used  by  MacCormack  and  Candler 
(Ref.  2)  to  calculate  the  nonequilibrium  chemistry  flow  about  a  two-dimensional  sphere- 
cone  at  Mach  18.  Test  conditions  for  the  calculation  were  as  follows: 


Mach  number: 

Freestream  Velocity: 
Static  Pressure: 
Static  Temperature: 


17.94 
5940  m/s 
79.779  N/ m2 
270. 65K 


(19488  ft/s) 
(1.67  lbf/ft2) 
(487°R) 


? 
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The  nose  radius  of  the  blunt  body  was  0.2  meters  (7.9  inches)  and  the  cone  section  had  a 
half-angle  of  10  degrees.  The  Reynolds  number  based  on  nose  radius  was  7.16  (104).  An 
isothermal  wall  temperature  of  1500K  (2700° R)  was  assumed,  and  five  gas  species  (N2, 
O2,  NO,  N,  O)  and  three  vibrational  energies  were  included.  Shock-fitting  was  used  to 
fit  the  strong  bow  shock.  The  results  presented  in  Figures  5  through  8  were  obtained  in 
approximately  150  implicit  steps  using  a  30x30  grid.  As  was  the  case  for  the  NS 
algorithm  test  cases  described  in  Section  2.2.6,  the  intent  here  was  to  assess  the  ability 
of  the  chemistry  algorithm  to  calculate  gross  flow  features  of  hypersonic  flows.  No 
attempt  was  made  to  validate  the  algorithm  to  confirm  the  accuracy  of  the  fine  flow 
details.  The  figures  presented  here  are  intended  only  to  demonstrate  global  calculation 
capability. 

Figure  5  shows  mass  concentrations  on  the  stagnation  streamline  for  the  four  minor 
species  (O2,  NO,  N,  and  O)  versus  axial  distance  along  the  missile,  while  Figures  6  (a-e) 
present  2D  contours  of  mass  concentration  for  each  of  the  five  gas  species. 
Temperature  ratios  on  the  stagnation  streamline  versus  axial  distance  along  the  model 
are  shown  in  Figure  7  for  the  translational-rotational  temperature  (T)  and  the  three 
vibrational  temperatures  for  N2,  NO,  and  O.  Figure  8a  presents  2D  contours  of  the 
translational-rotational  temperature,  while  Figure  8b  presents  2D  contours  for  the 
vibrational  temperature  of  N2. 
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Figure  6c.  Contours  of  Mass  Concentration  of  NO  (Percent) 
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Temperature  ratio 


(a)  Translational-Rotational  Temperature  Contours  (Kelvin) 
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(b)  Vibrational  Temperature  of  N2  Contours  (Kelvin) 


Figure  8.  Two-Dimensional  Temperature  Contours 


2.4  VECTORIZATION  AND  PARALLEL  PROCESSING 


2.4.1  Background 

The  goal  of  this  contract  is  to  accurately  and  efficiently  (i.e.,  one  hour  or  less) 
calculate  low  density  reacting  gas  flows  about  hypersonic  vehicles  150-200  nose  radii  in 
length.  Sections  2.2  and  2.3  have  described  in  detail  the  algorithms  which  will  be  used 
to  solve  the  coupled  3D  Navier-Stokes/chemistry  equation  set.  Adequate  flowfield 
resolution  in  boundary  layers,  near  shocks,  and  in  the  leeside  separation  regions  will 
require  a  large  number  of  grid  points  to  discretize  the  flowfield.  The  resulting 
computational  effort  will  probably  not  be  capable  of  meeting  the  stated  time  goal  on  a 
single  processor  machine.  Therefore,  it  is  required  that  parallel  processing  and 
vectorization  of  the  solution  algorithms  be  accomplished.  The  Phase  I  effort  in  this 
area  consisted  of  assessing  the  parallel  processing  and  vectorization  potential  for  the 
Navier-Stokes  and  chemistry  solution  algorithms.  A  two-fold  approach  was  chosen:  (1) 
assessment  of  current  state-of-the-art  in  vectorization  and  parallel  processing  hardware 
and  software,  and  (2)  demonstration  of  parallel  processing  and  vectorization  on  a 
prototype  Navier-Stokes  solution  algorithm.  These  topics  will  be  discussed  separately. 

2.4.2  Assessing  the  State-of-the-Art 

Recent  developments  in  the  supercomputing  environment  (in  the  past  10  years), 
have  rapidly  advanced  the  state-of-the-art.  Projections  for  the  next  five  years,  with 
respect  to  computing  speed  and  storage,  are  perhaps  even  more  ambitious.  Within  the 
past  two  years,  a  number  of  mini-supercomputers  have  been  developed,  in  addition  to 
continued  development  of  recognized  supercomputers  such  as  the  Cray  series  of 
computers. 

The  supercomputer  environment  can  be  described  by  two  factors:  (1)  computer 
hardware  (architecture,  processors,  storage),  and  (2)  software  systems.  These  two 
factors  are  not  mutually  exclusive,  as  the  software  systems  are  generally  highly 
dependent  upon  the  computer  system  architecture.  However,  the  discussion  here  will  be 


restricted  to  software  for  vectorization  and  parallel  processing  only,  and  the  two  topics 
will  be  described  independently. 

Hardware.  By  far  the  most  well-known  supercomputers  in  use  today  are  the  Cray 
series  of  computers.  The  Cray  machines  represent  a  number  of  different  computer 
architectures,  ranging  from  the  silicon  technology  of  the  Cray-1  to  the  gallium  arsenide 
technology  projected  for  the  Cray-3  (due  out  in  1992).  Each  architecture  strives  to 
minimize  clock  time  while  maximizing  computing  efficiency.  As  such,  there  are  trade¬ 
offs  to  be  made;  for  instance,  chaining  was  omitted  in  the  Cray-2  due  to  physical 
dimensional  limitations  of  the  capacitors  to  be  used  for  chaining.  Another  problem 
facing  both  the  Cray-XMP  and  Cray-2  is  memory  contention  between  multiple  CPUs. 
Both  of  these  problems  can  be  alleviated  by  means  of  local  memory  which  reduces  the 
number  of  fetches  and  stores  to  main  memory.  With  memory  contention  reduced,  Cray 
plans  to  implement  up  to  64  processors  on  the  Cray-3. 

The  success  achieved  by  the  Cray  series  of  computers  has  led  to  the  creation  of  a 
number  of  computer  companies  whose  goal  is  to  develop  mini-supercomputers.  These 
mini-supercomputers  seek  to  provide  computational  power  approaching  that  of  the  Cray 
computers  at  a  proportionately  lower  cost.  Many  of  these  machines  are  intended  to  be 
dedicated  single  user  or  multiple  (2  or  3)  user  systems,  with  the  hope  that  the 
throughput  from  the  mini-supercomputers  will  be  cost-effective  when  compared  to  the 
more  powerful  (but  frequently  overloaded)  Cray  machines.  Therefore,  a  realistic 
assessment  of  the  state-of-the-art  in  computer  hardware  should  include  a  study  of  the 
mini-supercomputers  available  today. 

The  mini-supercomputer  companies  have  shown  much  interest  in  the  prototype 
Navier-Stokes  code  upon  which  this  contract  effort  is  being  based.  To  date,  the 
prototype  code  has  served  as  a  benchmark  for  the  EPS,  ALLIANT,  and  INTEL  mini- 
supercomputers,  as  well  as  for  the  SCS-40  and  CONVEX  mini-supercomputers.  The 
latter  two  computers  are  currently  being  evaluated  in-house  by  the  Boeing  Aerospace 
Company  at  the  Kent  Space  Center.  Results  of  these  benchmark  tests  are  proprietary. 

The  architectures  of  the  FPS,  SCS-40,  and  CONVEX  computers  are  similar,  in  that 
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they  are  all  single  processor  machines.  The  ALLIANT  has  eight  processors  and  the  Intel 
computer  has  a  multiprocessor  hypercube  architecture  which  required  substantial  source 
code  modifications  before  the  benchmark  could  be  calculated.  All  of  the  mini¬ 
supercomputers  permitted  some  degree  of  vectorization.  However,  it  is  still  felt  that, 
for  the  purposes  of  the  current  contract,  the  target  machine  for  Phase  II  of  the  study 
should  be  a  Cray  series  computer.  The  Cray  system  computers  are  more  powerful  than 
the  mini-supercomputers,  and  have  a  more  established  softwa.e  system  for 
vectorization  and  parallel  processing  than  the  relatively  new  mini-supercomputers 
available  today. 

Software.  The  software  systems  on  the  aforementioned  mini-supercomputers  are 
currently  in  various  stages  of  development,  with  many  of  the  systems  (directly  or 
indirectly)  trying  to  mimic  the  software  systems  from  Cray  Research.  Therefore,  the 
discussion  here  will  focus  on  recent  developments  in  Cray  Research  software. 

Vectorization  still  represents  the  most  basic  method  for  improving  code 
performance,  and  is  executed  on  a  single  processor.  Vectorization  consists  of  arranging 
nested  DO  loops  such  that  the  longest  vector  occupies  the  innermost  loop.  The  only 
requirement  for  vectorization  is  that  all  DO  loop  index  dependencies  be  removed  on  this 
inner  (vector)  loop.  Vectorization  alone  can  result  in  significant  reduction  in  computer 
run  times. 

Additional  increases  in  computer  code  efficiency  can  be  realized  by  multi¬ 
processing  a  vectorized  code.  Cray  Research  currently  utilizes  two  techniques  to 
accomplish  parallel  processing-namely,  macrotasking  and  microtasking.  Macrotasking, 
which  consists  of  adding  compiler  directives  around  large  blocks  of  coding  that  may  be 
executed  simultaneously  on  multiple  processors,  has  existed  for  some  time  ar.d  is  well- 
documented.  Microtasking  is  a  more  recent  software  development,  and  consists  of 
adding  compiler  directives  around  small  blocks  of  coding  which  may  be  executed 
simultaneously  on  multiple  processors. 

Microtasking  is  basically  a  more  refined  version  of  macrotasking.  The  selection  of 
microtasked  blocks  is  more  complicated  than  the  selection  of  macrotasked  blocks; 
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however,  Cray  will  be  introducing  a  CFT77  compiler  (release  2.0  for  the  X-MP  and 
release  2.1  for  the  Cray-2)  which  will  automatically  microtask  (Ref.  17).  Additionally, 
Cray  will  be  introducing  a  new  compiler  directive  DO  GLOBAL  LONG  VECTOR,  which 
will  identify  long  vector  loops  for  special  treatment. 

2.4.3  Vectorization  and  Multiprocessing  Demonstration 

An  explicit  algorithm  version  of  the  prototype  three-dimensional  Navier-Stokes 
code  under  development  in  this  contract  was  vectorized,  and  a  test  case  was  run  on  a 
Cray  X-MP  using  a  single  processor.  A  flowtrace  of  the  solution  process  showed  that 
85%  of  the  CPU  time  was  spent  in  the  solution  module,  where  all  88  DO  loops 
vectorized. 

The  vectorized  code  was  then  multiprocessed  by  microtasking  the  explicit  solution 
module.  A  schematic  of  the  microtasking  operation  is  shown  in  Figure  9,  including 
compiler  directives  for  the  Cray  X-MP  COS  1.14  compiler.  The  microtasked  vision  of 
the  code  was  run  successfully  on  the  Cray  X-MP/24  (2  processors)  and  Cray  X-MP/48  (4 
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Figure  9.  Multiprocessing  the  Explicit  Algorithm  on  a  Cray  X-MP 
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processors).  Additionally,  the  test  case  was  run  on  1,  2,  4,  and  8  processors  on  the 
ALL1ANT  mini-supercomputer,  as  well  as  16  processors  on  the  Intel  Hype  s  be. 

Similar  vectorization  and  multiple  processing  operations  will  be  perfor  -d  on  the 
code  under  development  in  the  current  contract.  It  is  expected  that  with  vectorization 
and  multiple  processing,  the  contract  objectives  of  calculating  (in  less  than  1  hour)  the 
reacting  gas  flowfield  about  bodies  150-200  nose  radii  long  will  be  satisfied. 
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m.  AIR  CHEMISTRY  REACTION  MODEL  AND  WALL  CATALYSIS 


3.1  INTRODUCTION 

The  purpose  of  this  section  is  to  describe  the  work  performed  on  the  subtasks  for 
developing  an  air  chemistry  reaction  model  and  the  selection  of  reaction  rates  and 
equilibrium  constants.  Figure  10  summarizes  the  procedure  taken  in  this  phase  of  the 
contract. 


Figure  10.  Subtasks  for  Developing  a  Chemistry  Reaction  Model  and  the  Selection  of  Reaction  Rates 
and  Equilibrium  Constants 


The  sequence  of  steps  included  a  supplemental  literature  search,  construction  of  a 
detailed  air  chemistry  reaction  model,  a  sensitivity  analysis,  and  an  examination  of 
thermodynamic  and  transport  properties  methodology.  By  performing  these  tasks  a 
reduced  set  of  equations  containing  the  important  air  species  and  reaction  kinetics  was 
assembled.  Finally,  the  boundary  conditions  are  discussed. 
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3.2  DESCRIPTION  OF  GENERALIZED  EQUATIONS 

The  following  section  discusses  the  generalized  equations  that  can  be  used  to  model 
hypersonic  flows.  For  hypersonic  flow,  the  temperatures  produced  around  the  vehicle 
can  excite  rotation,  vibration,  chemical  ionization,  electron  excitation  modes,  and 
chemical  disassociation.  An  example  of  chemical  and  thermal  nonequilibrium  regimes 
about  the  stagnation  region  of  a  blunt  body  is  depicted  in  Figure  11. 


Figure  1 1.  Thermal  and  Chemical  Nonequilibrium  Regions  About  a  Stagnation  Point 

The  governing  equation  for  the  coupled  fluid  dynamic/chemistry  equation  set  can  be 
written  in  strong  conservation  form: 


au  aF  aG  aH 

_ c  +  _ c  ^  _ c  ^  _ c  _  ^ 

at  ax  ay  az 


(14) 


where  IFc  represents  the  conservative  solution  vector,  F*c»  G'c>  represent  the 
coupled  fluid  dynamic/chemistry  flux  vectors,  and  W  Represents  the  chemistry  source 
terms.  Primes  denote  cartesian  (nontransformed)  quantities. 

A  generalized  transport  equation  can  be  constructed  for  each  of  the  nonequilibrium 
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processes  which  comprise  the  chemical  source  term  in  Equation  (14)  (Ref.  18).  For  a 
given  chemistry  model  a  general  reaction  can  be  written  as 


N’s  k  V 

J3  t  r  ^ 

^  a  C  ^  ^  b  C  ,  r  = 

rs  s  rs  s 

s=  l  k  s  =  l 

br 


where  ars,  brs  are  the  stoichiometric  coefficients,  Cs  are  the  reactant  and  product 
species,  Ns,  Nr  are  the  number  of  species  and  reaction  paths,  and  kfr,  kbr  are  the 
forward  and  backward  reaction  rates,  respectively. 

The  law  of  mass  action  states  that  the  rate  of  change  of  concentration  of  species  s 
by  reaction  r  is  given  by 


11 1  c.  - k... 


The  rate  change  in  concentration  of  species  s  by  all  Nr  reactions  is  then  found  by 
summing  the  contributions  from  each  reaction. 


Cs  =  X  (C 


Finally,  the  production  rate  of  species  s  is  found  from 


=  C  m 
*  s  s 


The  forward  reaction  rates  are  computed  from  the  Arrhenius  Law: 


k  =  A  (T)  exp  1  -  E  /Ti 

f  r  r  r 


for  each  reaction  r.  The  coefficients,  Ar,  Nr  and  Er  are  experimentally  determined 
constants.  The  reverse  rates  can  be  found,  given  the  forward  rate  and  equilibrium 
constant  Keq  for  each  reaction,  as 


k  =  k.  /  K 

br  tr  eq 
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where  Keq  is  determined  from  the  Gibbs  free  energy  and  is  derived  in  Ref.  19. 

Thus,  it  can  be  seen  that  the  source  term  for  species  s  is  coupled  to  all  of  the 
reaction  paths  containing  specie  s.  The  reaction  equations  could  represent  O2 
disassociation  as 


o2^o  +  0 


(21) 


a  reaction  representing  a  diatomic  molecule  undergoing  a  transition  from  one 
vibrational  state  (VI)  to  another  vibrational  state  (V2)  as 


an  ionization  reaction  for  NO  as 


NO  ^  NO  *  -f-e~ 


(22) 


(23) 


or  an  electronic  excitation  reaction  for  O2  going  from  an  El  excited  mode  to  an  E2 
excited  mode 


O 


F. 

2 


1  __ 


(24) 


Reactions  of  these  types  can  exist  for  other  species.  In  fact,  the  number  of 
possible  reactions  needed  to  model  air  chemistry  is  very  large.  A  discussion  of  which 
reaction  equations  are  important  can  be  found  in  Section  3.3. 

Finally,  to  close  the  equation  set,  an  expression  for  the  total  internal  energy  is 
needed.  The  total  internal  energy  can  be  written  as 


e .  =  e.  +  e  +  e  ,  +  e  ,  +  e  .  +  e, 

1  trans  rot  vib  diss  elect  kin 


(25) 


where  etrans  represents  the  contribution  due  to  translation;  erot>  the  internal  energy 
due  to  rotation;  evib>  the  internal  energy  due  to  vibration;  eqiss,  the  internal  energy  due 
to  chemical  disassociation/association;  eeiect>  the  internal  energy  due  to  electronic 
states;  and  ekin  is  the  kinetic  energy  in  the  flow. 
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Each  of  these  energy  quantities  is  a  function  of  the  various  specie  concentrations 
and  may  not  be  in  equilibrium.  This  represents  the  general  partial  differential  equation 
set  describing  a  chemical  reaction  system.  Considerable  reduction  and  simplification  of 
the  equations  is  possible.  The  mechanism  for  reducing  the  complexity  of  the  problem  is 
discussed  in  Section  3.3. 

Typically,  if  all  of  the  internal  energy  modes  are  in  equilibrium,  then  the  total 
internal  energy  can  be  represented  by  a  polynomial  in  temperature.  Reference  20  gives 
an  example  of  how  this  is  implemented. 

Due  to  the  potential  importance  of  the  electron  distribution  in  the  flow  field, 
ionizational  effects  will  be  modeled.  Ionization  will  be  modeled  by  incorporating 
reaction  mechanisms  and  transport  equations  for  both  ions  and  electrons. 

One  simplification,  which  has  been  used  by  several  researchers  including  Vincenti 
(Ref.  21)  and  Wagener  (Ref.  22)  involves  combining  the  various  nonequilibrium  vibration 
modes  into  a  single  equation  for  the  nonequilibrium  vibrational  energy.  The 
simplification  results  by  assuming  that  a  diatomic  molecule  behaves  like  a  simple 
harmonic  oscillator.  With  this  assumption,  a  transport  equation  can  be  written  for  the 
nonequilibrium  vibrational  energy  as 


aF' 

c 

aG' 

c 

aH' 

r 

+  -  + 

at 

8x 

ay 

az 

Pev,b  -  pe, 


vi  b 


(26) 


where  Fc>  G'c>  and  EPc  are  defined  in  Equation  (14),  e*vib  is  the  equilibrium  vibrational 
energy  given  by 


N'vk 

e*  =  - 2- 

vib  0  /T 

e  v  -1 


(27) 


and  the  characteristic  vibrational  relaxation  time  scale  by 


l 


i  = 


(28) 


r 


f 


In  these  equations,  N  represents  the  total  number  of  oscillators,  8v  =  kpv/kB,  with  kp 
being  Planck's  constant,  v  the  frequency,  and  kB  Boltzmann's  constant.  Ki,o  is  a 
constant  dependent  on  the  particular  molecule. 

Another  physical  process  that  could  be  important  for  hypersonic  flow  analysis  is 
radiation.  The  incorporation  of  radiation  effects  involves  rewriting  the  energy  equation 
to  include  the  radiant  heat  flux  from  the  gas.  Applying  the  first  law  of  thermodynamics 
the  modified  energy  equation  can  be  written  as 


a(pe.)  aF'  aG'  aH'  aqR 

— — -  +  — -  +  — -  +  - -  =  -  i-  ,  (29) 

at  ax  ay  az  ax^ 


where  the  heat-flux  vector 
species  intensity  Iv  as 


can  be  expressed  in  terms  of  the  frequency-dependent 


f  00  f  4n 

qR  =  I  I  I  dQdv  ,  (30) 

J  Jo  Jo  VJ 

where  lj  is  the  unit  vector  specifying  the  direction  of  propagation  of  Iv  and  ll,  12,  13  are 
the  associated  direction  cosines.  lv  can  be  evaluated  using  the  equation  of  radiative 
transfer  in  the  quasi-equilibrium  form  as 


j 


(31) 


where  the  unsteady  term  has  been  dropped  due  to  the  fact  that  the  speed  of  light  is 
many  orders  of  magnitude  greater  than  typical  fluid  speeds  and  av  is  the  frequency- 
dependent  absorption  coefficient.  The  Planck  function  Bv  is  given  by 


BV(T)  = 


2  k  v3/c2 

p _ 

expl  vk  /k„ T1  - 1 


(32) 


where  kp  is  the  Planck  constant,  v  is  the  frequency,  c  is  the  speed  of  light,  and  kB  is  the 
Boltzmann  constant.  Thus,  the  energy  equation  is  described  by  a  complicated  integro- 
differential  relation.  This  equation  can  be  greatly  simplified  by  assuming  the  gas  to  be 
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thick  and  emission-dominated.  Vincenti  (Ref.  21)  outlines  how  these  simplifications  can 
be  implemented. 

In  addition,  both  translation  and  rotation  nonequilibrium  processes  tend  to  reach 
equilibrium  in  approximate  10  collisions  (Ref.  21).  Consideration  of  physical  events 
occurring  in  10  mean  free  paths  (10  collisions)  would  require  solving  the  noncontinuum 
equations.  Thus,  rotation  and  translation  can  always  be  modeled  as  being  in  equilibrium 
1  within  the  Navier-Stokes  approximation. 

The  various  conservation  equations  are  summarized  in  Figure  12.  These  equations 
can  be  solved  in  a  variety  of  w»vs  resulting  in  different  levels  of  nonequilibrium  physics. 
|  Several  combinations  of  these  equations  are  depicted  in  Figure  13. 


No. 

Conservation  equations 

No.  of  eq. 

What  quantities  are  obtained 

1 

Overall  mass  conservation  eq. 

1 

P 

n 

Overall  momentum  conservation  eq. 

3 

U,  V,  w 

Overall  energy  conservation  eq. 

1 

E  ■  "■►Ttrans 

■ 

Mass  conservation  eq.  for  each 
species 

Ns 

Pi 

5 

Vibrational  energy  conservation  eq. 
for  diatomic  molecules 

Nd 

evib,s - ►^vib.s 

6 

Electron  momentum  conservation  eq. 

3 

Electric  field ^ 

7 

Electron  energy  conservation  eq. 

1 

®elec - ►Te 

Mach  number  range  10-30 

>  Altitude  100,000  ft  -  non-continuum  limit  (-275,000ft) 

Figure  12.  Defining  Equations  for  Multi-Temperature  Hypersonic  Flows 


Model 

number 

Model 

type 

Temperatures 

modeled 

Conservation  equations 
(figure  12)  to  be  solved 

M-1 

NS 

Ttrans 

1.2,3 

M-2 

NS  +  CNE 

Ttrans  “  ^"rot  ■  "^vib- 

1.2, 3,  4 

M-3 

NS  +  CNE  +  TNE 

Ttrans  ”  "^rot-  "^"vib 

1.2,3, 4,5 

M-4 

NS  +  CNE  +  TNE 

Tfrans  =  Trot,  Tujb,  Tg^ 

1,2,3,  4,  5.6 

M-5 

NS  +  CNE  +  TNE 

Ttrans  •  Troj,  Tvib’  ^elec 

1,2, 3,  4, 5,  6, 7 

NS  —  Navier-Stokes 

CNE  —  Chemical  non-equilibrium 

TNE  —  Thermal  non-equilibrium 

Figure  13.  Equation  Modeling  Hierarchy 
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An  extensive  literature  search  was  performed  for  chemical  and  thermal 
nonequilibrium  measurements.  Several  sources  of  shock  tube  experiments  were  found 
(Refs.  23  and  24)  that  could  be  used  to  calibrate  nonequilibrium  air  chemistry 
mechanisms.  The  data  was  considered  sufficiently  accurate  to  compare  various 
chemical  nonequilibrium  models  (i.e.  one-temperature  models).  The  same  was  not  true, 
however,  for  thermal  nonequilibrium  (multi-temperature)  experimental  data  sets.  Very 
few  data  sets  were  found  which  contained  useful  measurements  of  thermal 
nonequilibrium  processes.  Based  on  the  results  of  the  literature  search,  it  was  decided 
that  only  one-temperature  air  chemistry  models  could  be  calibrated  to  engineering 
accuracy  and  would  be  considered  in  this  phase  of  the  contract.  Thus,  rotation, 
translation,  and  vibration  modes  would  all  be  assumed  to  be  in  equilibrium.  The 
conservation  equations  needed  to  model  this  type  of  flow  is  depicted  as  model  2  (M-2)  in 
Figure  13. 

Thermal  radiation  from  the  shock  front  can  be  significant  at  extremely  high  speeds 
(above  8-10  km/sec),  in  general,  thermal  radiation  from  the  shock  front  is  dependent  on 
the  body  bluntness  and  the  free-stream  velocity  and  density  of  the  gas.  The  radiation 
contribution  to  the  total  (convective  and  radiative)  stagnation  point  heat  transfer 
increases  with  velocity,  density,  and  body  bluntness  while  the  depth  of  the 
nonequilibrium  region  behind  the  bow  shock  varies  inversely  with  gas  density.  For  a 
blunt  (nose  radius  of  3  m.)  vehicle  traveling  at  7.6  km/sec,  the  convective  heating  will 
be  at  least  an  order  of  magnitude  larger  than  the  (equilibrium)  radiative  contribution 
above  an  altitude  of  46  km  (Ref.  25).  The  radiative  contribution  will  be  even  smaller  if 
this  vehicle  had  a  sharper  nose  geometry  and/or  traveled  at  higher  altitudes.  Hence, 
thermal  radiation  from  the  shock  front  can  be  assumed  to  be  small  relative  to 
convective  heat  transport  mechanisms  and  was  not  modeled. 

3.3  CHEMICAL  REACTION  MODEL  DEVELOPMENT 

A  summary  of  the  chemical  reaction  model  development  is  depicted  in  Figure  14. 

Details  of  all  stages  of  the  development  process  are  documented  below. 


Figure  14.  Summary  of  Air  Chemistry  Model  Development  Process 

3.3.1  Literature  Search 

An  extensive  literature  search  for  air  ionization  models,  specie  thermodynamic 
properties  (enthalpy,  specific  heat  capacity,  entropy,  and  Gibb's  free  energy)  and  specie 
experimental  measurements  was  performed  during  Phase  I.  Many  air  ionization  models 
were  identified,  with  each  model  utilizing  different  chemical  species,  reaction  paths, 
rate  coefficients,  and  application  features.  For  instance,  a  model  may  depend  on  one 
(translation)  or  two  (translation  and  vibration)  temperatures,  may  have  up  to  15 
chemical  species,  may  use  an  Arrhenius  or  non-Arrhenius  reaction  rate,  and  may  have  as 
many  as  50-60  reaction  paths.  These  variations  make  a  thorough  analysis  extremely 
tedious  and  difficult.  To  insure  a  feasible  analysis,  the  following  assumptions  and 
restrictions  were  employed: 

a.  only  Arrhenius  reaction  rates  were  considered; 

b.  chemical  reaction  rates  were  assumed  to  depend  only  on  the  translational 
temperature  (T),  i.e.,  one  temperature  reaction  rates; 
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c.  chemical  reaction  rates  were  assumed  to  be  valid  up  to  17000  K; 

d.  argon  has  a  negligible  effect  in  high  temperature  air  reactions. 

Of  the  many  models  examined  only  five  air  ionization  models  utilized  an  Arrhenius 
reaction  model  of  the  form  Ar(T)^rexp  (-Er/T)  where  Ar,Nr,Er  are  rate  coefficients. 
These  models  were  developed  by  K.  Wray  (hereafter  abbreviated  by  W  ),  C.  Park  and 
G.  Menees  (PM),  S.  Kang  and  M.  Dunn  (KD),  D.  Bittker  (Bit),  and  M.  Bortner  (Bort) 
(Refs.  26-30).  The  W  model  is  well  known  ana  represents  one  of  the  simplest  (18 
reaction  paths,  7  species)  ionization  models  conceived.  The  PM  model  was  derived  to 
investigate  meteoroid  wakes  and  both  the  KD  and  Bort  models  were  employed  in 
flowfield  analyses  about  reentry  vehicles.  Finally,  the  Bit  model  was  employed  in  Ref. 
29  to  analyze  shock-tube  air  reactions.  The  reaction  paths  and  rate  coefficients  for 
each  model  are  depicted  in  Tables  1  through  5  and  summarized  in  Table  6.  Note  that 
the  W  model  depicted  in  Table  1  has  been  modified  from  its  original  form. 
Modifications  included  the  elimination  of  N2  +  02  =  2  NO,  based  on  Camac  and 
Feinberg's  recommendation  (Ref.  23)  and  the  assumption  that  argon  has  a  negligible 
effect  in  high  temperature  air  reactions. 


Table  1.  Wray  Air  Model 


Reaction 

M 

Nr 

Er 

O2  +  M  — *  20  +  M 

N,  NO 

n2 

°2 

O 

3.62  x  10'  9 

4.8  x  10  20 

1.9  x  1021 

6.4  x  10  23 

-1.0 

-1.5 

-1.5 

-2.0 

59500 

59500 

59500 

59500 

N2  +  M->2N  +  M 

O,  02,  NO 

n2 

N 

1 .9  x  1o]Z 

4.8  x  1 0 1 ' 

2? 

4.1  x10^ 

-0.5 

-0.5 

-1.5 

113000 

113000 

113000 

NO  +  M— >N  +  0  +  M 

N2,  02 

NO,  0,  N 

3.9  x  10  20 
7.9x1021 

-1.5 

-1.5 

75500 

75500 

NO  +  O  ->  02  +  N 

3.2x109 

1.0 

19700 

N2  +  O  ->  NO  +  N 

7.0x  1013 

0.0 

38000 

N  +  O  ->  NO+  +  e' 

6.5X1011 

0.0 

31900 

M  represents  the  collision  specie  in  a  given  reaction 
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Table  2.  Park  and  Menees  Air  Model 


Reaction 

M 

N, 

Er 

O2  +  M  — *  0  +  0  +  M 

N 

8.25  x  id9 

-1.0 

59,500 

0 

8.25  x  101 9 

-1.0 

59.500 

n2 

2.75  xIO19 

-1.0 

59,500 

02 

2.75  xIO19 

-1.0 

59,500 

NO 

2.75  xIO19 

-1.0 

59,500 

N2  +  M-*N  +  N  +  M 

N 

1.11  xIO22 

-1.6 

113,200 

0 

1.11  xIO22 

-1.6 

113,200 

N2 

3.7  xIO21 

-1.6 

113,200 

02 

3.7  xIO2' 

-1.6 

113.200 

NO 

3.7x1021 

-1.6 

113,200 

NO  +  M  _♦  N  +  0  +  M 

N 

4  6  xt  O1 7 

-0.5 

75,500 

0 

4.6x1017 

-0.5 

75.500 

N2 

2.3  xIO17 

-0.5 

75.500 

o2 

2.3x1017 

-0.5 

75.500 

NO 

2.3  xIO17 

■0.5 

75,500 

NO  +  O-.N  +  02 

2.16x10® 

1.29 

19,220 

O+N2  -+  N  +  NO 

O.ISxlO13 

0.10 

37,700 

N  +  0-»NO+  +  e- 

1.53  xIO10 

0.37 

32.000 

O  +  6  — ¥  C*  t  «  +  0 

1.95  xIO34 

-3.78 

158,500 

N  +  e  — *  N*-  +  e'+  e 

1.25  xIO35 

-3  82 

168,600 

0  +  0-»02*+e‘ 

3.85  xIO10 

0  49 

80,600 

0  +  0 2+-*02  +  0+  - 

6.85  xIO13 

-0  52 

18,600 

N2  +  N-i — *  N  +  Ng * 

S  85  xIO*2 

-0.18 

12,100 

N  +  N  — ►  N2+  +  e" 

1.79  xIO10 

0.77 

67,500 

0+  N0+-+N0  +  O+ 

2.75  xIO13 

0.01 

51,000 

N2  +  Of  -+ 0+  N^ 

6.33  xIO13 

-0.21 

22.200 

N  +  NO+-*NO+N+ 

2.21  xIO15 

-0.02 

61,100 

O2  +  NK>  -*  NO  +  02+ 

1.03  xIO16 

-0  17 

32,400 

N0++N-+N  +  N2* 

1.70  xIO1 3 

0.40 

35.500 

M  represents  the  collision  specie  in  a  given  reaction 


5' 


Table  3.  Dunn  and  Kang  Air  Model 


Reaction 

M 

Nr 

Er 

O2  +  M  — *  20  +  M 

N,  NO 

3.6  x  1018 

-1.0 

59400 

N2  +  M  -+2N  +  M 

O,  NO,  02 

1.9x  1017 

-0.5 

113000 

NO  +  M  -» N  +  0+  M 

O5,  N? 

3.9x  1020 

-1.5 

75500 

O+NO  ->N  +  02 

3.2x  10  9 

1.0 

19700 

o  +  n2  ->n  +  no 

7. Ox  1013 

0.0 

38000 

N  +  N2  -»2N  +  N 

4.085x  10  22 

-1.5 

113000 

0  +  N  -*NO"'  +  e' 

1.4x  10  6 

1.5 

31900 

O  +  e"  ->  0+  +  2  e~ 

3.6  x  10  31 

-2.91 

158000 

N  +  e”  — >N+  +  2e‘ 

1.1  x  10  32 

-3.14 

169000 

0+0  -»02++  e" 

1.6x  1017 

-.98 

80800 

o+o2+  -+o2+o+ 

2.92  x  1018 

-1.11 

28000 

n2  +  n+  -+n  +  n2+ 

2.02x  1011 

0.81 

13000 

N  +  N  — » N2++  e~ 

1.4x  1013 

0.0 

67800 

0+  NO+  -» NO  +  0+ 

3.63x  1015 

-.6 

50800 

n2  +  o+  -+o  +  n2+ 

3.4  x  10  19 

-2.0 

23000 

N  +  NO+  — >  NO  +  N+ 

I.Ox  1019 

-.93 

61000 

o2  +  no+  -+no  +  o2+ 

1.8x  1015 

0.17 

33000 

0+  NO+  -»02  +  N+ 

1.34  x  10 13 

0.31 

77270 

02  +  O  — *  20  +  O 

9.0  x  101 9 

-1.0 

59500 

02  +  02  — » 20  +  02 

3.24x1019 

-1.0 

59500 

02  +  N2  — >20  +  N2 

7.2x1018 

-1.0 

59500 

N2  +  N2  — »  2N  +  N2 

4.7  x  1017 

-0.5 

113000 

NO+M  -+N  +  O  +  M 

0,  N,  NO 

7.8  x  10  20 

-1.5 

755000 

02  +  N2  -+  NO  +  NO+  +  e‘ 

1.38  x  10  20 

-1.84 

141000 

NO+N2  -+NO+  e”  +  N2 

2.2  x  1 0 13 

-0.35 

108000 

NO  +  O2  — ♦  NO*  +  -f  O2 

: 

i 

8.8  x  1015 

-0.35 

108000 

M  represents  the  collision  specie  in  a  given  reaction 
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Table  4.  Bittker  Air  Model 


Reaction 

M 

Ar 

Nr 

ErR 

n  +  o2  ->no  +  o 

6.4  x  109 

1.0 

6250 

0  +  N2  -»NO  +  N 

1.8x  1014 

0.0 

76250 

N  +  O  +  M  — >  NO  +  M 

All  but  N20 

6.4  x  1016 

-0.5 

0 

n  +  o  +  n2o  ->no  +  n2o 

1 .44  x  1 0 1 7 

-0.5 

0 

O  +  O  +  M  ->  02  +  M 

All 

5.7  x  10 19 

0 

-1788 

N  +  N  +  M  ->  N2  +  M 

All 

2.8  x  1017 

-0.75 

0 

NO  +  O  +  M  ->  N02  +  M 

All  but  N2 

5.62X1015 

0 

-1160 

no  +  o  +  n2  -»no2  +  n2 

8.71  xIO15 

0 

-1160 

n2o  +M  ->N2  +  0  +  M 

All 

1.42x  1014 

0 

51280 

o  +  n2o  — »n2  +  o2 

6.23  XIO13 

0 

24520 

NO+  +  e*  ->N  +  0 

1.45x  1021 

-1.5 

0 

0+  +  e'  +  M  -»0  +  M 

All  but  N,  02,  NO,  O 

2.0  x  10  26 

-2.5 

0 

0+  +  e'+  N  -»0+  N 

6.0  x  10  24 

-2.5 

0 

O  •+•  a  +  C>2  —*  O  -*■  (j'2 

9.0  x  10  29 

-2.5 

0 

0+  +  e'  +  NO  — *  O  +  NO 

1. Ox  10  28 

-2.5 

0 

O  +  +  e~  +  O  — >  O  +  O 

6.0  x  10  24 

-2.5 

0 

02  +  e'  +  M  — » 02”+  M 

All  but  N2 

1. 52  x  10  21 

-1.0 

1190 

02  +  e~  +  N2  — >  02~  +  N2 

3.04  xIO1 6 

-1.0 

1190 

02  +  0  — »  02  +  O 

6.0  x  1012 

0 

0 

M  represents  the  collision  specie  in  a  given  reaction 
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Table  5.  Borlner  Air  Model 


Reaction 

Ar 

Nr 

Er 

02  +  M  — >  O  +  O  +  M 

m-o2 

1,37  x  10  '5 

-.83 

59400 

O 

1 .50  x  10  "4 

-1 

59400 

n2 

2  2  x  10-3 

-1.7 

59400 

N,  NO 

2.4  x  10  ‘3 

-1.8 

59400 

n2+  m->  n  +  n  +  m 

M-  N2 

8.3  x  10  '6 

-.75 

113200 

N 

5. Ox  10  -2 

-1.5 

113200 

02,  0,  NO 

4.8  x  10  '6 
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113200 

NO  +  M->  N  +  O  +  M 

M-  NO,  O,  N 

1.32  x  10'2 
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02.  ^2 

1.5  x  10  ‘9 

0 
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O3  +  M  — >  O  +  O2  +  M 

M  =  N2 
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0 
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19100 

O2  +  O2  — *  0  +  O3 

2.9  x  10  ‘11 

0 
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N  +  0-»  NO+  +  e‘ 

8.6x  10  _13 

0 

31900 

N2  +  M  ->  N2+  +  e'  +  M 

3.0  x  10 ‘20 

M  =  N2 

2.0 

181000 

e~ 

5.3x  10  -11 

0.5 

181000 

O2  +  M  — >  O2 4  +  e  +  M 

I.Ox  10  -20 

M  -  N2 

1.5 

140000 

o2 

2.0  x  10  '20 

1.5 

1 40000 

e” 

2.2x  10  -11 

0.5 

1 40000 

02  +  e”  +  M  — >  02  +  M 

4.2  x  10  '21 

M  =  02 
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600 

N2 

1  x  10  '21 

0 

0 

02  +  O3  —>  02  +  O3 

3  x  1 0  "10 

T/300 

0 

02"  +  O  — >  O3  +  e" 

2.5  x  10  '10 

T/300 

0 

O3  +  e"-)0't  02 

4x  10'11 

T/300 

0 

X++ V-+X  +  Y 

6  x  1 0'5 

-1 

0 

X+Y"  +  M->X  +  Y  +  M 

3  x  1 0 "7 9 

-2.5 

0 

N2+  +  O2  — >  N2  +  C>2+ 

1  x  lO-10 

T/300 

0 

N2+  +  NO-+N2  +  NO+ 

5x10-10 

T/300 

0 

n2+  +  o-+no+  +  n 

2.5x  10  '10 

T/300 

0 

o2+  +  no^o2  +  no+ 

O 

O 

X 

CO 

T/300 

0 

M  represents  the  collision  specie  in  a  given  reaction 
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Table  6.  Summary  of  Species,  Paths,  and  Rate  Constants  for  Various  Air  Models 


Reference 

Species 

Number  of 
reaction  paths 

Forward 
reaction  rate 

K.  L.  Wray 

O,  02,  N,  N2,  NO,  NO+,  e" 

18 

Ar  (Tpr  expf-Er/T) 

D.  A.  Bittker 

O,  0+,  O',  02‘,  02,  N,  N2,  NO, 

NO+,  N20,  N02,  e' 

89 

Ar(T)Nrexp(-Er/T) 

M.  H.  Bortner 

O,  O".  O2,  O2”,  02+,  03, 03" 

N,  N+,  N",  N2,  N2+,  NO,  NO+ 

>38 

Ar  (T)Nr  exp  (-Er/T) 

C.  Park  & 

G.  P.  Menees 

O,  0+,  02,  02++,  N,  N+,  N2,  N2+, 
NO,  NO+,  e' 

29 

Ar  (T)Nr  exp(-Er/T) 

M.  G.  Dunn  & 

S.  W.  Kang 

O,  0+,  02, 02+,  N,  N+,  N2,  N2+ 

NO,  NO+,  e” 

32 

Ar  (T  )Nr  exp  (-E/T) 

where  Ar,  Nr, 
T 

Er  *  reaction  rate  constants 
-  temperature 

Specie  thermodynamic  properties,  specifically  entropy,  Gibbs  free  energy,  enthalpy, 
and  specific  heat  capacity  are  important  components  to  model  development.  These 
properties  must  be  accurately  known  up  to  very  high  temperatures  (15000  K-25000  K). 
Two  sources  for  thermodynamic  properties  were  found.  D.  Esch  et  al,  (Ref.  31)  derived 
polynomial  relationships  for  the  thermodynamic  properties  of  99  elemental  and 
molecular  enpciec  More  than  one  half  of  all  the  species  were  modeled  up  to  15000  K. 
Similarly,  J.  Shinn  (Ref.  32)  obtained  polynomial  relationships  for  thermodynamic 
properties  of  11  air  ionization  species  (0,02>N,N2,N0,0+,02+,N+,N2+,N0+,e')  up  to 
30000  K.  Both  Esch  and  Shinn  utilized  the  following  relationships  for  the  specific  heat 
capacity  (Cps),  enthalpy  (hs),  entropy  (Ss),  and  Gibbs  free  energy  (fs)  of  specie  s: 


Cps/R  =  Zis  +  Z2sT  +  Z3sT2  +  Z4sT3  +  Z5sT4  ,  (33) 

hs/RT  =  Zls  +  Z2sT/2  +  Z3sT2/3  +  Z4sT3/4+ Z5sT4/5  +  Z«S,T,  (34) 

S«/R  =  Zl3ln(T)  +  Z2sT  +  Z3sT2/2  +  Z4sT3/3+ Z5sT4/4  +  Z7s,  (35) 

f</RT  =  Zls[l.-ln(T)]  -  Z2sT/2  -  Z3sT2/6-  Z4sT3/12  -  Z5sT4/20  +  Z6s/T  -  Z7s.  (36) 


A  careful  examination  of  Shinn's  curve-fit  coefficients  revealed  that  inaccurate 
results  are  obtained  for  N2  near  13500  K  -  14990  K.  Specifically,  negative  N2  heat 
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capacities  were  calculated  using  coefficients  Zis  through  Z5S.  Coefficients  Zis-Z5S  are 
also  used  to  estimate  entropy,  enthalpy,  and  free  energy,  hence,  these  properties  are 
also  unreliable.  No  anomaly  was  identified  for  the  remaining  10  species.  Because  of  the 
discrepancies  in  Shinn's  data,  Esch's  thermodynamic  data  was  used  in  the  model 
development. 

Near  the  end  of  Phase  I  a  third  source  of  thermodynamic  properties  was  obtained 
from  R.  Gupta  of  NASA  Langley  Research  Center.  Gupta's  curve-fit  coefficients 
models  the  previously  mentioned  11  air  species  up  to  25000  K  in  the  same  form  as 
Equations  (33)  through  (36).  Unfortunately,  it  was  not  made  available  in  time  to  use  in 
the  development  of  the  chemistry  reaction  model,  however,  his  results  will  be  examined 
in  Phase  II. 

Not  all  species  of  interest  have  been  modeled  to  high  temperatures  (15000  K-25000 
K).  Molecular  species  N2O  and  NO2  used  by  Bittker  can  be  computed  only  up  to  6000  K 
while  species  O3  and  O3-  considered  by  Bortner  have  not  been  modeled  at  all.  Due  to 
the  lack  of  adequate  thermodynamic  data  for  these  species  both  Bittker  and  Bortner 
models  were  not  analyzed. 

Finally,  experimental  measurements  are  needed  to  validate  the  air  models  at  high 
speeds.  An  extensive  literature  search  identified  only  two  useful  sources  of  specie 
measurements.  Camac  and  Feinberg  (Ref.  23)  recorded  the  formation  of  nitric  oxide 
(NO)  in  shock  heated  air  for  Mach  Numbers  6-11.  Peak  NO  concentrations  and  time-to- 
peak  periods  were  measured  in  this  Mach  number  range.  For  Mach  numbers  above  14, 
electron  density  measurements  behind  a  shock  front  were  investigated  by  Lin  (Refs.  24, 
33).  Complete  electron  density  histories  were  measured  for  Mach  Numbers  14-20. 

3.3.2  Nitric  Oxide  Calculations 

In  this  part  of  the  analysis  the  accuracy  of  the  W,  PM,  and  KD  models  were 
assessed  in  a  comparison  of  nitric  oxide  calculations  with  shock  tube  measurements  by 
Camac  and  Feinberg  (Ref.  23).  Peak  NO  predictions  and  measurements  are  depicted  in 
Figure  15.  The  W,  PM,  and  KD  models  yielded  results  within  the  measurement  band. 
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Figure  15.  Comparison  of  KD,  W ,  and  PM  Model  Calculations  With  Peak  Nitric  Oxide  Measurements 
of  Camac  and  Fein  berg  ( Reference  23) 

The  predicted  time-to-peak  periods  are  revealed  in  Figure  16.  All  models  predicted 
much  higher  peak  periods  at  a  shock  velocity  of  2.5  mm/miero-sec.  The  differences 
between  measurements  and  predictions  may  be  due  to  Camac  and  Feinberg's 
interpretation  of  the  peak  time.  Infrared  radiation  histories  (from  which  NO 
concentrations  were  obtained)  at  this  shock  velocity  have  longer  molecular  relaxation 
periods  resulting  in  a  "flat"  peak  history  that  can  be  difficult  to  interpret. 
Nevertheless,  all  models  predicted  identical  results  at  this  velocity.  Calculations  at 
higher  velocities  indicated  that  the  W  and  PM  models  yielded  good  agreement  with 
measurements.  Examination  of  Figures  15  and  16  show  that  the  W  and  PM  models  are 
suitable  for  high  temperature  NO  calculations. 


3.3.3  Electron  Density  Calculations 

W,  PM,  and  the  KD  electron  density  histories  downstream  of  a  shock  wave  are 
compared  with  Lin  et.  al.  (Ref.  24)  shock  tube  measurements  in  Figures  17  through  19 
for  Mach  15.6,  18.5,  and  20.2.  All  results  were  computed  with  the  original  Esch 
thermodynamic  properties.  Figure  17  shows  all  models  underpredicted  peak  and 

M) 


Figure  17.  W,  PM,  and  KD  Electron  Density  Calculations  (particles /cm° I  Versus  Distance,  (cm) 
Downstream  of  a  Shock  Front  for  Mach  15.6  Using  Esch  Thermodynamic  Properties 


Figure  19b.  W  and  PM  Electron  Density  Calculations  ( part  ides /cm^i  Versus  Distance  (cm)  Downstream 
of  a  Shock  Front  for  Mach  20.2  Using  Esch  Thermodynamic  Properties 

equilibrium  electron  measurements  by  approximately  40%  at  Mach  15.6  and  identify  the 

importance  of  obtaining  accurate  equilibrium  predictions.  At  Mach  18.5  (Fig.  18)  the 

KD  model  overpredicted  peak  electron  measurements  by  approximately  4  times  while 

both  W  and  PM  results  underestimated  measurements.  The  differences  between  models 

is  more  evident  at  Mach  20.2  depicted  in  Figure  19  (a  and  b).  The  KD  peak  electron 

densities  are  now  25  times  larger  than  measurements,  the  PM  results  are  50%  larger, 

and  again  the  W  model  underpredicted  measurements  by  70%.  The  W  model  yielded 

results  that  were  consistently  lower  than  the  measurements  at  all  shock  velocities.  All 

models  predicted  earlier  time-to-peak  periods  at  Mach  18.5  and  20.2. 

The  results  presented  in  Figures  17  through  19  illustrate  the  importance  of 
obtaining  accurate  equilibrium  electron  densities.  The  equilibrium  concentration  is 
governed  by  the  equilibrium  constant  (Keq)  which  can  be  defined  in  terms  of  entropy  and 
enthalpy, 
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where  Ns  =  the  number  of  species,  ars,  brg  =  stoichiometric  coefficients  of  the  sth 
specie  in  the  rth  reaction.  The  latter  expression  was  obtained  from  polynomial 
relationships  for  entropy  and  enthalpy  mentioned  earlier.  Note  that  Keq  can  be 
modified  by  varying  the  "Z7S"  coefficient. 

Improved  results  were  obtained  using  the  (modified)  Esch  thermodynamic  properties 
when  the  Z7S  coefficient  for  the  electron  density  was  increased  by  10%  from  its  original 
value.  These  results  are  shown  in  Figures  20  through  22.  All  models  showed  significant 
improvements  at  Mach  15.6  (Fig.  20).  A  10%  increase  in  (electron)  Keq  resulted  in  a 
45%  increase  in  peak  and  equilibrium  concentrations.  At  Mach  18.5  (Fig.  21)  the  KD 
model  continued  to  overpredict  measurements  at  all  locations  downstream  of  the  shock 
front.  The  PM  model  adequately  predicted  equilibrium  concentrations  but  overshoots 
peak  experimental  data.  Satisfactory  results  were  obtained  with  the  W  model  at  all 
stations.  Electron  histories  at  Mach  20.2  are  depicted  in  Figure  22  (a  and  b).  The  PM 
and  KD  models  yielded  results  much  higher  than  measurements.  The  W  model 
accurately  reproduced  the  peak  concentration  and  the  time-to-peak  period  but 
underpredicted  the  near-equilibrium  value. 
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ire  20.  I V,  PM,  and  KD  Electron  Density  Calculations  (particles /cm^)  Versus  Distance  (cm) 

Downstream  of  a  Shock  Front  for  Mach  15.6  Using  Modified  Esch  Thermodynamic  Properties 


RUN  DATA©  W 
.RUN  DATA©  PH 
RUN  DATA©  KD 

'run  data©  Meas 


n,  run,  ana  is  u  tiectron  uensity  calculations  I particles  /cm''/  Versus  Distance  (cm) 
Downstream  of  a  Shock  Front  for  Mach  1&5  Using  Modified  Esch  Thermodynamic  Properties 
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Figure  22a.  W,  PM,  and  KD  Electron  Density  Calculations  (particles /cm^)  Versus  Distance  (cm) 

Downstream  of  a  Shock  Front  for  Mach  20.2  Using  Modified  Esch  Thermodynamic  Properties 
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Figure  22b.  W  and  PM  Electron  Density  Calculations  (particles /cm^)  Versus  Distance  (cm)  Downstream 
of  a  Shock  Front  for  Mach  20.2  Using  Modified  Esch  Thermodynamic  Properties 
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3.3.4  Sensitivity  Analysis 

W  and  PM  analytical  results  depicted  in  Figures  20  through  22  are  rearranged  in 
Figures  23  and  24  for  Mach  15.6,  18.5,  and  20.2.  Although  the  W  and  PM  models 
continue  to  yield  unsatisfactory  results  they  represent  a  feasible  beginning  for  the 
sensitivity  analysis. 

The  W  model  satisfactorily  predicted  peak  concentrations  but  underpredicted  near¬ 
equilibrium  concentrations  above  Mach  15.6  (Fig.  23).  The  W  model  has  only  one 
reaction  path  directly  responsible  for  electron  production  and  several  non-ionized  paths 
that  can  have  an  indirect  effect  (Table  1).  Perturbation  of  the  W  model  rate 
coefficients  (Ar,Nr,Er)  by  one  or  two  orders  of  magnitude  failed  to  yield  satisfactory 
peak  and  near-equilibrium  values. 

The  PM  model,  on  the  other  hand,  yielded  satisfactory  equilibrium  concentrations 
but  predicted  higher  peak  densities  above  15.6  (Fig.  24).  The  PM  model  (Table  2)  has 
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Figure  24.  PM  Electron  Density  Calculations  ( particles /cmr)  Versus  Distance  (cm)  Downstream  of  a 

Shock  Front  for  Mach  15.6.  18.5,  and  20.2  Using  Modified  Esch  Thermodynamic  Properties 

electron  production,  hence,  it  is  more  likely  to  yield  improved  results  in  a  sensitivity 
analysis.  Perturbation  of  several  reaction  rate  coefficients  yielded  significant 
reductions  in  peak  values  while  having  little  effect  on  equilibrium  calculations.  For 
instance,  a  5%  decrease  in  the  "Nr"  coefficient  for  the  0  +e-=  2e_  +  0+  and  N  +e~=  2e~ 

+  N+  reaction  paths  resulted  in  a  45%  reduction  in  peak  electron  concentrations  at  Mach 
20.2.  Much  smaller  peak  reductions  were  observed  at  lower  Mach  numbers.  Larger 
reductions  in  the  "Nr"  coefficient  were  ineffective.  Perturbation  of  the  rate 
coefficients  for  all  other  paths  had  little  effect  on  the  results.  Finally,  further 
reduction  of  peak  values  were  obtained  by  deleting  paths  21  through  29  in  Table  2. 
These  paths  generated  excessive  quantities  of  electrons  at  the  higher  Mach  numbers  and 
had  little  effect  on  the  non-ionized  species  (0,02,N,N2,N0)  .  This  version  of  the  PM 
model  (Table  7)  yielded  the  best  overall  agreement  with  electron  density  measurements 
up  to  Mach  20.2  (Fig.  25)  using  the  modified  Esch  thermodynamic  properties. 
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Table  7.  Modified  Park  and  Menees  Air  Model 


Modified 

Park  and  Menees 
Air  Model 


M  represents  the  collision  specie  in  a  given  reaction 

*  Originally  -  3.78 

#  Originally  -  3.82 

3.3.5  Transport  Properties 

The  transport  of  specie  momentum,  thermal  energy,  and  mass  is  characterized  by 
the  coefficient  of  viscosity  (ps),  thermal  conductivity  (\s),  and  (concentration)  diffusion 
(Dst)>  respectively.  Transport  properties  can  have  a  considerable  impact  on 
nonequilibrium  flow  calculations.  For  instance,  viscosity  affects  the  growth  of  the 
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Figure  25.  Modified  PM  Electron  Density  Calculations  t, particles /err, Versus  Distance  (cm)  Downstream 
of  Shock  Front  for  Mach  15.6,  18.5,  and  20.2  Using  Modified  Esch  Thermodynamic  Properties 


boundary  layer,  hence,  influences  the  distance  (and  time)  particles  and  energy  must 
travel  to  reach  the  wall  surface.  Conversely,  thermal  conductivity  regulates  the 
transfer  of  heat,  thus,  influences  particle  motion  (i.e.  specie  concentration).  Finally,  it 
will  be  shown  that  moderately  accurate  prediction  methods  are  typically  costly  and 
lengthy. 

All  chemical  nonequilibrium  methods  require  the  computation  of  each  transport 
property  for  individual  species  in  the  reactive  flowfield.  A  formula  is  then  used  to 
obtain  the  mixture  value  of  the  property.  This  procedure  is  depicted  in  Figure  26.  Also 
shown  is  the  functional  dependence  of  each  property  on  other  parameters.  Ii  is  worthy 
to  note  that  all  properties  depend  on  the  collision  cross-section  02fi<lml  which  is  defined 
as  the  probable  collision  area  between  species  s  and  t.  The  collision  cross-section  is 
dependent  on  specie  mass,  velocity,  temperature,  collision  deflection  angle,  and  an 
assumed  potential  force  function.  Details  can  be  found  in  Ref.  34. 
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where  p  - 

A  « 
0  - 
s  - 

mix  m 

M  . 
T  - 
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Cp  - 

P  - 

c  - 


viscosity 
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chemical  specie  s 
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collision  cross-section  area 
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Figure  26.  Typical  Transport  Properties  Methodology 

Specifically,  first-approximations  to  the  specie  transport  properties  can  be 
represented  as  follows  (Ref.  34): 


P 


S 


V  M  T 

266.93  X  10  ~7 - - - 

a2  Q<2-2’*  (T) 
s  s 


(38) 


A 

S 


=  1989.1  X  10“7 


\/T/M 

s 


o2ft 

s 


(2,2)* 

s 


(T) 


(39) 


pD  ,  =  2.628  X  10~3 

St 


\/T3(M  +  M  )/2  M  M 

_ S _ t _ s  t 

02qIUI*(T) 
s  St 


(40) 


where  o2Q,lmr  is  the  collision  cross-section  normalized  to  its  rigid-sphere  value. 
Equations  (38)  to  (40)  are  based  on  kinetic  theory  of  dilute  gases  and  complete 
derivations  are  provided  in  Ref.  34. 

Simplifications  to  Equations  (38)  to  (40)  have  been  developed  by  Tong  (Ref.  35). 
The  transport  properties  were  curve-fitted  with  an  equation  of  the  form 
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(41) 


Q  =  exp  [a.  (In  1  8T)2  +  Ps  (In  1 ,8T)  +  ys  ] 

where  Q  represents  ps,  \s,  pDst  and  as,  Ps»  Ys  are  the  curve-fit  coefficients.  The 
coefficients  for  each  specie  property  are  listed  in  Tables  8  through  17.  Equation  (41)  is 
accurate  to  within  0.5%  of  the  source  data  which  was  generated  from  a  computer 
program  developed  by  Hatch  (Ref.  36)  and  Pindroh  (Ref.  37).  The  source  data  is  valid  up 
to  20,000K  and  this  method  will  most  likely  reduce  computation  time  and  storage 
requirements.  Note  that  Tong's  method  does  include  collisions  with  non-ionized  and 
ionized  species  but  has  neglected  coulomb  effects  (i.e.  electron  collisions).  The  effect 
of  neglecting  coulomb  collisions  can  be  evaluated  during  the  Phase  II  validation  test 
trials. 

Once  the  specie  transport  properties  are  known  a  mixture  value  for  each  property 
may  be  obtained  by  one  of  several  formulas.  An  excellent  review  of  mixture 
relationships  for  thermal  conductivity  and  viscosity  are  discussed  in  Refs.  38  and  39, 
respectively.  In  general,  nonequilibrium  mixture  formulas  involve  the  summation  of  the 
product  (and/or  quotient)  of  the  specie  concentration  and  its  transport  property  for  each 
specie  in  the  flow.  One  of  the  simplest  and  straightforward  method  to  compute  the 
thermal  conductivity  is  a  linear  mixing  formula  (Refs.  38,  40), 
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Table  8.  Specie  Dynamic  Viscosity 


Specie 

a 

P 

y 

0 

.03044 

.2242 

-11.04 

02 

-.01718 

.9345 

-13.69 

NO 

-.01578 

.9147 

-13.69 

N 

.03811 

.09672 

-10.49 

NO+ 

-.0153 

4.856 

-41.10 

n2 

.0110 

-16.  14 

50.39 

N+ 

-.1530 

4.856 

-41.48 

o+ 

-.1530 

4.856 

-41.41 

In  ns  =  a  (In  1 ,8T)2  +  (5  (In  1 ,8T)  +  y 


Table  9.  Specie  Thermal  Conductivity 


Specie 

a 

P 

y 

O 

.03045 

.2242 

-11.80 

o2 

-.01718 

.9346 

-15.15 

NO 

-.Cl  578 

.9147 

15.08 

N 

.03811 

.09674 

-11.12 

NO* 

-.1530 

4.856 

-42.49 

n2 

1.097 

-16.  14 

49.07 

N+ 

-.1530 

4.856 

-42.11 

o+ 

-.1530 

4.856 

-42.17 

In  =  a  (In  1.8T)2  +  p  (In  1.8T)  +  y, 

X  (  _££! - - -  ).TfK) 

'  cm-sec-K ' 


L 
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Table  10.  Binary  Diffusion  Coefficient  of  Specie 
0  Into  Specie  s 


Specie 

a 

0 

O 

.03344 

1.18750 

-9.55591 

o2 

.04139 

1 .07857 

-965397 

NO 

.04060 

1.09031 

-9.70358 

N 

.01866 

1.50710 

-11.19786 

NO+ 

.04060 

1.09031 

-9.70358 

n2 

.05342 

.88732 

-8.94594 

N  + 

.06302 

.75995 

-9.16222 

o+ 

.00169 

1.55148 

-11.69451 

In  (pDos)  =  P  18T)2  +  0  (In  1.8T)  +  y 

_crn2\  ,  p(3tm),  T  (°K) 
sec  ' 

Table  12.  Binary  Diffusion  Coefficient  of  Specie 
NO  Into  Specie  s 


Specie 

a 

0 

y 

0 

.040598 

1.09031 

-9.70358 

o2 

-.000376 

1.69266 

-12  20127 

NO 

.0017725 

1.66381 

-12  10346 

N 

-.0097108 

1.82599 

-12  33594 

NO  + 

.0032396 

1.56533 

-12  42772 

n2 

.022085 

1.35095 

-10.94333 

N  + 

-.0097108 

1.82598 

•12.33594 

o+ 

,040598 

1.C9031 

Table  1 1.  Binary  Diffusion  Coefficient  of  Specie 
O 2  Into  Specie  s 


Specie 

at 

P 

y 

O 

.04139 

1 .07857 

-9.65397 

02 

-  005643 

1.77269 

-12.49255 

NO 

-.000376 

1 .69266 

-12.20127 

N 

-.009046 

1.81422 

-12.27983 

NO+ 

-.000376 

1.69266 

-12.20127 

n2 

-.002632 

1.64120 

-11.99645 

N  + 

-.0090458 

1.81422 

-12.27983 

o+ 

.04139 

1.07857 

-9.65397 

In  (pDo2s  )  =  cx  (In  1.8T)2  +  0  (In  1.8T)  +  y 

Q02s(_2£2)  ,  p(atm),  T  (°K) 


Table  13.  Binary  Diffusion  Coefficient  of  Specie 
N  Into  Specie  s 


Specie 

tx 

p 

y 

O 

01 6664 

1.50710 

-11.19786 

°2 

-.0090458 

1.81422 

-12.27983 

NO 

-.0097108 

1.82599 

-12.33594 

N 

.028369 

1.30193 

-10.05754 

NO  + 

-.0097108 

1.82599 

-12.33594 

N2 

.075705 

.563218 

-7.77918 

N  + 

.001265 

1.55013 

-11.64697 

O1' 

-.063017 

.75995 

-9.16222 

In  (pDf^Qg)  «  a  (In  1.8T)2  +  0  (In  1.8T)  +  y 

dNOs(  cm2  ^  ■  P(atm)'  T  (  K) 

'  sec  I 

Table  14.  Binary  Diffusion  Coefficient  of  Specie 
NO+!nto  Specie  s 


Specie 

a 

P 

y 

O 

.04060 

1.09031 

-9.70358 

02 

-.000376 

1.69266 

-12.20127 

NO 

.0032394 

1.56533 

-12.42772 

N 

-.0097108 

1.82599 

-12.33594 

NO+ 

.0032394 

1.56533 

-12.42772 

n2 

.022085 

1.35095 

-10.94333 

N+ 

-.0097108 

1.82599 

-12.33594 

o+ 

.040598 

1.09031 

-9.70358 

In  (pDiMa)  «  u  (In  1.8T)2  +  0  (In  1.8T)  +  y 

DNs(  cm2  \  ,  p(atm),  T  (=K) 

'  sec  I 

Table  15.  Binary  Diffusion  Coefficient  of  Specie 
N 2  Into  Specie  s 


Specie 

a 

P 

y 

O 

.053422 

.887319 

-8.94594 

°2 

.002632 

1.64120 

-11.99645 

NO 

.022085 

1.35095 

-10.94333 

N 

.075705 

.563218 

-7.77918 

NO+ 

.022085 

1.35095 

-10.94333 

n2 

1.12800 

-15.  58970 

52.74193 

N  + 

7.57048 

.563218 

-7.77918 

0+ 

.053422 

.887319 

-8.94594 

In  (pC^0+s)  =  a  (In  1  8T)2  +  0  (In  1 .8T)  +  y 

DNO+sf  cm2  )  ■  p(a,m)’  T  (°K) 

'  sec  • 


In  (pDfsj^g)  =  a  (In  1 ,8T)2  +  0  (In  1 .8  T)  +  y 

Dn  (  cm2  \  .  p(atm).  T  (CK) 

“  '  sec  I 


Table  16.  Binary  Diffusion  Coefficient  of  Specie 
N+  Into  Specie  s 


Specie 

a 

f» 

Y 

0 

.063017 

.759948 

-9.162217 

°2 

-.0090458 

1.81422 

-12.27983 

NO 

-.0097108 

1 .825985 

-12.33594 

N 

.001265 

1.550133 

-11.64697 

NO  + 

-.0097108 

1.825985 

-12.33594 

N2 

.075705 

.563218 

-7.77918 

N  + 

.0012645 

1.55013 

-11.64697 

0+ 

.063017 

.759948 

-9.162217 

n  :  pCfj  *s  )  =  a  (In  1 ,8T)2  +  p  (In  1 ,8T)  +  y 


DN*s  •  P(alm).  TTK) 

'  sec  ' 


Table  17.  Binary  Diffusion  Coefficient  of  Specie 
0+  Into  Specie  s 


Specie 

a 

P 

Y 

O 

.001690 

1.55148 

-11.69451 

o2 

.041390 

1.07857 

-9  15397 

NO 

.040598 

1.09031 

-9.70358 

N 

.0630167 

.759948 

-9.16222 

NO  + 

.040598 

1.09031 

-9.70358 

n2 

.053422 

.887319 

-8.94594 

N  + 

.0630167 

.759948 

-9.16222 

o+ 

.0016898 

1.55148 

-1 1.69451 

In  (pDQ+g)  =  u  (In  1.8T)2  +  [5  (In  1.8T)  +  y 


Dq+o  (  cm2  \  ,  P(atm).  T  (°K) 
'  sec  I 


Ns  C  (42) 

\  =  M  X  —  .\  , 

mix  mix  —  s 


where  Cs 


mass  fraction  of  specie  s, 


\s  =  thermal  conductivity  of  specie  s, 
Ms  =  molecular  weight  of  specie  s, 
Mmix  =  mixture  molecular  weight. 


Results  obtained  from  Equation  (42)  are  typically  larger  than  experimental 
measurements  (Refs.  38,  41-43).  On  the  other  hand-  the  reciprocal  mixing  rule, 


,  C  (43) 

-  -  M  N  — —  , 

\  —  M  \ 

mix  >  -  1  -x  s 


yield  results  smaller  than  measured  values  (Ref.  43).  Subsequently,  Burgoyne  and 
Weinberg  (Refs.  38,  42)  recommended  combining  Equations  (42)  and  (43)  to  obtain 


\ 

nn\ 


\1 


N  . 

V 


C  t 

s  s 

M 

S 


1 


s 

\1  v 

in;  v  — 

s  -  I 


c 

s 

M  \ 

>  s 


(44) 
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Equation  (44)  has  been  validated  for  rare  gases  (Ref.  44),  Ar-He  systems  (Ref.  45),  and 
polyatomic  gases  (Ref.  46)  with  moderate  accuracy  (less  than  20%  deviations).  The 
simplicity  of  this  expression  makes  it  attractive  for  numerical  programming. 

For  diatomic  and  polyatomic  molecules  As  must  include  effects  due  to  internal 
degrees  of  freedom  not  present  in  monatomic  species.  A  modified  Eucken  expression 
for  the  thermal  conductivity  is  suggested  (Ref.  47), 

Cps  ms  (45) 

A  =  .352  -  +  .12  A  , 

s.E  l  R  s’ 

where  Cps  =  specie  specific  heat  capacity, 

R  =  universal  gas  constant. 

The  specie  thermal  conductivity  can  be  shown  to  be  directly  proportional  to  the 
specie  viscosity  (based  on  first-order  approximations,  Ref.  34),  hence,  a  relationship 
analogous  to  Equation  (44)  can  be  obtained  for  the  dynamic  viscosity, 


The  accuracy  of  Equation  (46)  is  expected  to  be  similar  to  its  thermal  conductivity 
counterpart. 

Finally,  the  reciprocal  of  the  mixture  diffusion  coefficient  has  been  found  to  be 
approximately  equal  to  the  mass  average  of  the  reciprocal  of  the  binary  coefficient 
(Refs.  48,  49),  that  is, 

Ns 

v  C 

D  =  - IfLifJ-* _ 

mlx-1  Ns  c  (47) 

M  V  _JL_ 
mnt  “  MD 

S  =  1 .  S  X  t  S  St 

Individual  species  formulas  38-40  ana  mixture  Equations  (44),  (46),  and  (47)  are  the 
recommended  relationships  for  determir;ng  the  transport  properties.  These  equations 
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may  be  utilized  several  million  times  in  a  complete  three-dimensional,  Navier-Stokes 
calculation  about  a  simple  geometry.  The  transport  property  process,  therefore, 
requires  considerable  computational  effort  and  cost-effectiveness/accuracy  will  be 
carefully  assessed  during  Phase  11.  If  necessary  a  scheme  to  discriminate^  employ  these 
relationships  at  selected  conditions  (temperature  ranges,  locations)  may  be  utilized  to 
reduce  computational  effort.  These  equations  also  preclude  the  effects  of  chemical 
reactions  (dissociations,  recombinations)  on  the  transport  properties.  Chemical 
reactions  have  a  significant  effect  on  the  total  thermal  conductivity  and  a  moderate 
influence  on  concentration  diffusion  at  verv  high  temperatures  (Ref.  37).  To  include 
these  effects  would  require  additional  lengthy  calculations  and  increase  the  total 
computation  time.  Nevertheless,  the  effectiveness  of  the  above  method  to  predict  heat 
transfer  and  electron  density  histories  can  only  be  assessed  during  Phase  II  test  trials. 

3.4  WALL  CATALYSIS 

In  this  task  the  appropriate  boundary  conditions  were  developed  to  account  for  wall 
catalytic  effects.  Details  of  this  procedure  will  depend  on  the  numerical  algorithm 
selected.  In  the  nonequilibrium  flow  over  a  body  at  hypersonic  speeds,  the  body  surface 
may  act  as  a  catalyst  for  the  recombination  of  atoms  and  ions,  hence  increasing  the 
heat  transfer  at  the  surface.  Reentry  heating  data  derived  from  STS-2,  STS-3,  and  STS- 
5  missions  clearly  showed  the  significance  of  nonequilibrium  gas  chemistry  on 
aerodynamic  heating  in  a  high-velocity,  low-density  flight  regime.  As  altitude  and/or 
velocity  is  increased  above  STS  reentry  levels,  nonequilibrium  effects  will  become  even 
more  pronounced.  Consequently,  having  the  capability  to  treat  chemical  reactions, 
including  surface  kinetics,  is  essential  in  developing  a  flowfield  code  for  treating  high¬ 
speed,  low-density  flows. 

The  chemistry  conditions  considered  at  the  wall  include  non-catalytic,  finite  rate, 
and  fully-catalytic  processes.  Each  condition  is  summarized  below  with  reference  to 
Figure  27: 
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Figure  27  Specie  Catalysis  at  the  Wall  Surface 


Non-catalvtic 

ri  -  p 

diss  "  '-'diss 

c  -  c 
^ rec  '-'rec 


Finite  rate 


1  -  A  P 

diss  ~  "diss^diss’ 

0 

<  ^diss 

1  -  Q 

rec  ~  Drec^rec’ 

0 

*  Brec 

Fullv-catalvtic 
^'diss  =  ® 

C'rec  =  ^reelstoichiometric 


where  C 

^diss’  Brec 
diss 

ree 


specie  concentration 

constants  indicating  the  extent  of  finite  rate  wall  catalysis 
dissociated  species 
recombined  species 


At  least  two  different  types  of  thermal  boundary  conditions  will  be  modeled.  The 
first  case  involves  specifying  a  constant,  uniform  wall  temperature  and  the  second  case 
assumes  the  wall  to  be  adiabatic.  In  both  cases  the  degree  of  wall  catalicity  is  variable. 
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IV.  LEESIDE  EFFECTS  AND  TURBULENCE  MODELS 


4.1  LEESIDE  EFFECTS 

4.1.1  Statement  of  Work 

The  objective  of  the  leeside  model  task  was  to  develop  an  improved  flowfield  model 
for  leeward  surfaces  of  hypersonic  vehicles  at  angle-of-attack.  An  accurate 
determination  of  the  aerotherma.  environment  in  these  regions  is  needed  to  ensure  an 
adequate  thermal  protection  system/structural  design,  and  to  avoid  the  cost  and 
performance  penalties  associated  with  over-design. 

4.1.2  Background 

When  vehicle  angle-of-attack  is  sufficient  to  separate  the  flow  over  leeside 
surfaces,  vortices  are  formed  that  can  strongly  influence  heating  patterns.  These 
leeside  flows  are  highly  sensitive  to  vehicle  geometry,  angle-of-attack,  and  Re* t.oids 
number,  and  have  proven  very  difficult  to  analyze. 

To  minimize  empiricism  in  the  present  analysis  (and  thus  increase  reliability  and 
generality),  the  Phase  I  leeside  flow  analysis  attempted  to  calculate  leeside  flows 
directly,  without  adding  any  empiricism  over  and  above  that  introduced  by  the 
turbulence  models  to  be  discussed  in  Section  4.2.  Minimizing  empiricism  was  considered 
important  since  the  sensitivity  exhibited  by  leeside  flows  studied  to  date  suggests  that 
direct  extrapolations  to  flight  conditions  may  not  be  very  reliable.  An  additional 
advantage  of  the  current  (direct)  approach  is  that  it  facilitates  the  inclusion  of  chemical 
reaction  effects  in  the  leeside  vortices. 

4.1.3  Technical  Development 

The  standard  Baldwin-Lomax  turbulence  model  (Ref.  50)  has  been  widely  applied  to 
computing  supersonic  turbulent  flowfields  around  conical  bodies  at  low  and  moderate 
angles  of  attack.  Good  agreement  between  computed  results  and  experimental  data  can 
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be  obtained  except  in  the  leeside  region  where  vortex  flow  separation  occurs.  The 


reasons  for  this  discrepancy  are  twofold:  improper  determination  of  the  turbulence 
length  scale  and  poor  grid  resolution  incapable  of  handling  the  leeside  vortex  flow 
structure. 

The  major  difficulty  encountered  in  the  leeside  computation  is  the  determination  of 
the  proper  length  scale  upon  which  to  base  the  change  in  the  eddy  viscosity  calculation 
from  the  formula  valid  in  the  inner  region  to  that  in  the  outer.  In  the  original  Baldwin- 
Lomax  model  this  length  scale,  ymax>  was  determined  as  the  y  value  for  which  F(y) 
below  was  a  maximum: 

F(y)  =  wy  [1-exp  ( —  y +/261] ,  (48) 

where  <j  is  the  vorticity,  y  is  the  distance  normal  to  the  wall,  and  y+  is  the  law-of-the- 
wall  coordinate.  In  the  leeside  region  vortex  flow  separation  causes  F(y)  to  have  two  or 
more  maxima.  The  original  Baldwin-Lomax  formulation  would  choose  the  second 
maximum,  yielding  an  improperly  high  value  of  outer  eddy  viscosity.  Degani  and  Sehiff 
(Ref.  51)  suggested  a  modification  to  the  Baldwin-Lomax  model,  choosing  invariably  the 
first  maximum  of  F(y)  in  the  leeside  region  (Fig.  28).  Using  grids  with  adequate 
resouticn  for  the  leeside  vortices,  they  applied  the  modified  Baldwin-Lomax  model  to 
supersonic  flows  around  various  ogive-cylinder  bodies  and  cones  at  high  angle-of-attack. 
The  computed  results  thus  obtained  were  in  good  agreement  with  the  experiment  data 


Figure  28.  Behavior  of  Leeside  Function  F  (Y)  at  High  Angle-of-Attack 
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even  in  the  leeside  region.  Figures  29,  30  show  some  typical  results  obtained  by 
Rainbird  (Ref.  52)  using  the  modified  Baldwin-Lomax  formulation. 


Figure  29.  Normalized  Pitot-Pressure  Contours  in  Flowfield  for  a  Yawed  5-deg  Cone 
(Mx  =  1.8,  a  =  12. 5- deg,  Rex  =  28.9  (l(fi}) 


Figure  30.  Computed  Crossflow  Plane  Velocity  Vectors  on  a  5-deg  Cone 
(Mx  =  1.8,  a=  12. 5- deg,  Re x  =  28.9  ( 1(fi 1/ 

It  is  concluded  that  using  the  Baldwin-Lomax  model  with  Degani-Shiff 
modifications  is  adequate  to  compute  the  flowfield  on  the  leeside  region  of  a  body  at 
various  angles  of  attack  and  mach  numbers,  and  should  be  sufficient  to  satisfy  the 
leeside  calculation  requirements  of  the  present  contract. 

Part  of  the  work  performed  for  this  task  included  a  literature  review  to  identify 
data  sets  that  could  be  used  to  evaluate  the  leeside  calculation  capability  of  the  code  to 
be  developed  during  Phase  II.  Six  candidate  cases  have  been  identified,  and  are  listed  in 
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Table  18  below.  These  eases  span  a  wide  range  of  Mach  number/angle-of-attack 
combinations,  and  should  provide  an  adequate  supply  of  test  results  to  compare  against. 
These  cases  will  be  supplemented  during  Phase  II  if  necessary,  and  if  new  test  data  is 
identified. 


Table  18.  Candidate  Leeside  Evaluation  Cases 


Configuration 

Flight  condition 

Reference 

o,  deg 

Mach 

S’  sharp  cone 

12  5 

1  8.4  25 

52 

5  6°  sharp  cone 

2  to  18 

14  2 

53 

9° sharp  cone 

0  to  20 

70 

54 

15°  sphere  cone 

■5  to  15 

5  25,  7  40,  10  6 

55 

12  8/7° biconic 

-10  to  40 

60 

56 

13*17°  bicomc 

Oto  10 

100 

57 

4.2  TURBULENCE  MODELS 

4.2.1  Statement  of  Work 

The  high-altitude,  low-density  (and  thus  low  Reynolds  number)  flows  to  be 
considered  in  the  current  contract  will  be  primarily  laminar  flows.  In  the  late  phase  of 
re-entry,  however,  transition  to  turbulence  wi’l  occur.  Therefore,  a  contract 
requirement  was  to  include  at  least  one  turbulence  model.  This  model  was,  in  principal, 
separate  from  the  leeside  turbulence  model. 

4.2.2  Background 

The  turbulence  expected  for  the  case  of  interest  in  this  contract  will  include  shear 
layer  type  turbulence  on  the  windward  side  of  the  re-entry  vehicle  as  well  as  vortex 
separated  flows  on  the  leeward  side  of  the  vehicle.  Windward  turbulence  of  the  shear 
layer  type  is  expected  to  occur  as  a  shock/boundary  layer  interaction  at  a  deflected 
control  surface,  for  example.  This  shear  layer  type  turbulence  is  relatively  simple  to 
model,  even  with  basic  algebraic  turbulence  models  such  as  the  Baldwin-Lomax  (Ref.  50) 
model.  For  the  Phase  I  work  described  here,  both  algebraic  and  two-equation  turbulence 
models  were  studied.  Leeside  models  were  discussed  in  Section  4.1. 
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4.2.3  Technical  Development 

Extensive  literature  review  indicated  that  the  standard  Baldwin-Lomax  eddy 
viscosity  turbulence  model  (Ref.  50)  has  been  successfully  used  in  wide  ranges  of 
applications  in  computational  fluid  dynamics.  This  algebraic  turbulence  model  is 
relatively  simple  to  implement  and  gives  accurate  results  in  transonic  and  supersonic 
flows  without  separation.  Its  accuracy  deteriorates  somewhat  in  flows  with  separation 
but  nonetheless  it  still  gives  acceptable  results  in  many  engineering  applications.  As  the 
first  step  of  this  task  the  standard  Baldwin-Lomax  turbulence  model  was  installed  into 
the  prototype  three-dimensional  explicit  Navier-Stokes  code.  The  implementation  was 
validated  under  independent  research  and  development  funding  using  two-dimensional 
viscous  transonic  flows  around  airfoils.  Excellent  agreement  between  the  computed 
results  and  wind  tunnel  data  were  obtained  (Ref.  58). 

In  order  to  study  the  option  of  a  more  sophisticated  turbulence  model  for  flows  in 
the  hypersonic  Mach  number  range,  a  literature  review  was  carried  out  with  emphasis  on 
two-equation  models.  In  general,  two-equation  turbulence  models  show  a  slight 
improvement  over  the  standard  Baldwin-Lomax  model  for  computing  the  flowfields  at 
the  flight  conditions  prescribed  in  this  contract.  The  slight  improvement  in  solution 
thus  gained,  however,  may  not  be  very  cost-effective  due  to  the  nature  of  the  expected 
flowfields  and  the  severe  computational  cost  penalty  associated  with  two-equation 
formulation. 

One  candidate  turbulence  model  for  possible  inclusion  in  the  code  is  the  Jones- 
Launder  model.  This  model  was  introduced  in  1972  (Ref.  59)  and  has  proved  very 
popular  among  two-equation  models.  The  Jones-Launder  model  adds  one  partial 
differential  equation  for  the  turbulent  kinetic  energy  and  one  for  the  energy  dissipation 
rate.  A  comparison  of  the  Jones-Launder  model  with  other  turbulence  models  for  many 
of  the  expected  flows  mentioned  above  is  reported  by  Viegas  and  Horst  man  (Ref.  60).  In 
this  paper,  the  Jones-Launder  model  was  compared  with  the  Wilcox-Rubesin  two- 
equation  model  (Ref.  61),  the  Viegas  one-equation  model  (Ref.  62),  and  the  Baldwin- 
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Lomax  (Ref.  50)  model  for  a  variety  of  flow  problems.  As  expected,  the  two-equation 
models  generally  yielded  better  results. 

The  recent  approach  of  adding  wall  functions  to  the  turbulence  models  has 
circumvented  the  need  for  clustering  many  grid  points  near  solid  surfaces.  The  reduced 
grid  clustering  permits  longer  time  steps  to  be  taken  and  improves  computational 
efficiency.  A  typical  wall  function  was  suggested  by  Viegas  and  Rubesen  (Ref.  63). 
Results  of  calculations  using  this  wall  function  with  the  Jones-Launder  model  compared 
to  those  using  the  Baldwin-Lomax  model  for  a  compression  ramp  test  case  are  reported 
by  Knight,  et.  al.  (Ref.  64).  The  combination  of  the  two-equation  model  together  with 
the  wall-function  boundary  conditions  yields  slightly  better  results  in  general.  The 
computational  time  required,  however,  is  higher  and  the  complexity  of  the  code  will  be 
greatly  increased.  Therefore,  it  is  not  recommended  that  a  two-equation  turbulence 
model  be  added  to  the  Navier-Stokes  code  under  development. 


V.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  Phase  I  effort  of  the  current  contract  was  intended  to  identify  components  of  a 
combined  fluid  dynamic/chemistry  solution  algorithm  for  calculating  low  density  real 
gas  flows  about  hypersonic  vehicles.  It  is  believed  that,  by  utilizing  results  of  related 
Air  Force  funded  research  and  Boeing  project  and  research  work,  the  Phase  1  study  was 
more  thorough  than  originally  planned.  In  particular,  results  from  Air  Force  contracts 
AFOSR-85-0372  and  F33315-86-C-3015  supporting  research  work  at  Stanford  University 
under  Professor  R.  W.  MacCormack,  and  Air  Force  contract  F04611-86-C-0015  awarded 
to  Boeing  Aerospace  Company  for  the  development  of  an  advanced  Navier-Stokes  rocket 
base  flow  calculation  capability,  have  proved  most  useful  for  the  current  contract 
effort.  The  combined  results  of  these  supporting  contracts  and  the  development  work 
performed  in  Phase  I  of  this  contract  will  provide  a  strong  basis  for  the  Phase  II  program 
development  effort. 

At  the  start  of  this  contract  it  was  believed  that,  in  the  present  computer  hardware 
environment,  software  portability  would  play  a  major  role  in  the  long-range  success  of 
the  code  to  be  developed.  Experience  gained  during  the  Phase  I  effort,  in  particular 
during  the  evaluation  of  vectorization  and  parallel  processing  requirements,  has  served 
to  reinforce  the  importance  of  portable  software  systems. 

Sections  Il-IV  have  discussed  in  detail  the  conclusions  derived  from  the  Phase  I 
technical  effort  in  the  areas  of  solution  algorithm  and  chemistry  model  development, 
turbulence  modeling,  leeside  modeling  and  wall  catalysis.  The  recommendations  derived 
from  the  Phase  I  conclusions  for  application  to  the  Phase  II  program  development  effort 
are  summarized  below: 

•  Target  Computers:  Cray  X-MP/Y-MP  and  Cray  2/3 

•  Multiprocessing:  Cray  microtasking 
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•  Algorithm: 


Fully  coupled  chemistry 
Implicit  formulation  with  flux  splitting 
Gauss-Seidel  line  relaxation  solution  procedure 

•  Chemistry:  Modified  Park  and  Menees  air  chemistry  model  with  3 

species  and  20  reaction  paths 

Modified  Esch  thermodynamic  properties  up  to  15,000  K 
Tong  simplification  of  Pindroh  model  for  transport 
properties 

•  Turbulence  Model:  Modified  Baldwin-Lomax  model 

Basing  the  Phase  II  program  development  effort  upon  the  above  recommendations  will 
result  in  computer  program  capable  of  accurately  and  efficiently  predicting  trends  in 
the  heating  rates,  pressures,  forces,  and  moments  on  complex  shapes  in  low-density 
reacting  gas  hypersonic  flows. 
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ABBREVIATIONS 


ADI 

alternating  direction  implicit 

AFB 

Air  Force  Base 

BAC 

Boeing  Aerospace  Company 

Bit 

Bittker  air  chemistry  model 

Bort 

Bortner  air  chemistry  model 

CFD 

computational  fluid  dynamics 

FNS 

full  Navier-Stokes 

GS 

Gauss-Seidel 

KD 

Kang-Dunn  air  chemistry  model 

NS 

Navier-Stokes 

PM 

Park-Menees  air  chemistry  model 

PNS 

parabolized  Navier-Stokes 

RFP 

request  for  proposal 

RNS 

reduced  Navier-Stokes 

TLNS 

thin  layer  Navier-Stokes 

TVD 

total  variation  diminishing 

W 

Wray  air  chemistry  model 

2D 

two-dimensional 

3D 

three-dimensional 
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NOMENCLATURE 


Math  Symbols 


a 

= 

partial  derivative 

D+ 

- 

forward  difference 

D- 

- 

backward  difference 

D± 

- 

alternating  forward  and  backward  differences 

exp 

- 

exponential 

In 

= 

natural  logarithm 

At 

= 

time  increment 

GO 

= 

infinity 

£ 

= 

summation 

8Un+1 

= 

implicit  change  in  solution  vector 

AUn 

= 

explicit  change  in  solution  vector 

— 

= 

forward-backward  chemical  reaction 

-> 

= 

forward  chemical  reaction 

— 

= 

vector 

/ 

= 

integral 

V 

= 

square  root 

Ns 

n 

s-\ 

= 

product  of  terms  with  index  from  s  =  1  to  s  = 

Units 

m/s 

= 

meters  per  second 

ft/s 

= 

feet  per  second 

N/m2 

Newtons  per  square  meter 

lbf/ft2 

= 

pound-force  per  square  foot 

K 

Kelvin 
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°R 

- 

degrees  Rankine 

cm 

= 

centimeters 

cal 

= 

calories 

sec 

s 

seconds 

atm 

= 

atmospheres 

km 

= 

kilometers 

m 

= 

meters 

Subscripts 

i.  j,  k 

= 

coordinate  index  discretizations 

r 

= 

r-th  reaction 

s,  t 

= 

species  designators 

inv 

= 

inviscid 

vis 

= 

viscous 

max 

= 

maximum 

Superscripts 

i 

+ 


n 

T 

A 


refers  to  physical  (nontransformed)  flux  vectors 
positive  eigenvalue  contributions  (except  for  y+  in  Equation  48) 
negative  eigenvalue  contributions  (except  for  e~  electron 
designation) 

time  level  discretization 

vector  transpose 

refers  to  diffusion  velocities 
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Greek  Symbols 


av 


as,  Ps,  Yg 


i,,  n,< 

Sx>  Sy»  Sz 
lx,  1y,  Hz 
£x»  4y»  £z 
°2st^st^*  m)* 


©V 

aa,  ab,  ac 

\ 

P 

v 

P 

t; 

« 

<j) 


absorption  coefficient  (Equation  31) 

curve-fitted  coefficients  for  the  specie  transport  properties 
(Equation  41) 

general  curvilinear  coordinates 

metrics  of  the  general  curvilinear  transformation 

metrics  of  the  general  curvilinear  transformation 

metrics  of  the  general  curvilinear  trpnsformation 

collision  cross-sectional  area  between  species  s  and  t,  normalized 

to  its  rigid  sphere  value  (Equations  38-40) 

characteristic  vibration  temperature  (Equation  28) 

diagonal  matrices  of  the  eigenvalues  of  A,  B,  C,  respectively 

thermal  conductivity 

dynamic  viscosity 

bulk  viscosity  =  -Jp 

frequency  (Equation  27) 

total  density 

shear  stress 

vorticity 

specie  net  production  rate 


English  Symbols 
arsi  brs 

*0,  «1»  *2,  *3, 
*4,  *5,  *6 

A,  B,  C 


stoichiometric  coefficients  of  the  reactants  and  products, 
respectively,  for  specie  s  in  the  r-th  reaction  path  (Equation  15) 
coefficient  matrices  for  the  linearized  3D  Navier-Stokes  finite 
difference  algorithm 

Jacobians  of  the  flux  vectors  F,  G,  H,  respectively 
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Ar>  Np, 

Ep 

= 

Arrhenius  rate  for  r-th  reaction  path  (Equation  19) 

Bv 

= 

Planck's  function  (Equation  32) 

c 

= 

speed  of  light 

C 

= 

specie  concentration 

CP 

= 

specific  heat  capacity  at  constant  pressure 

Dst 

= 

(concentration)  diffusion  coefficient  of  specie  s  into  specie  t 

e~ 

- 

electron 

e 

- 

total  energy  per  unit  volume  =  p[ej  +  i(u2+v2+w2)] 

ei 

- 

total  internal  energy  per  unit  mass  =  ediss+eelect+skin+epot+ 

etrans+evib 

ediss 

- 

internal  energy  per  unit  mass  due  to  dissociation 

Select 

- 

internal  energy  per  unit  mass  due  to  electrons 

®kin 

- 

internal  energy  per  unit  mass  due  to  kinetic  energy 

Srot 

- 

internal  energy  per  unit  mass  due  to  rotation 

strans 

- 

internal  energy  per  unit  mass  due  to  translation 

Svib 

- 

internal  energy  per  unit  mass  due  to  vibration 

e*vib 

- 

equilibrium  vibrational  internal  energy  per  unit  mass 

f 

- 

Gibb's  free  energy  (Equation  36) 

F(y) 

= 

function  in  the  Baldwin-Lomax  turbulence  model 

F,  G,  H 

= 

x,  y,  z  momentum  flux  vectors  for  the  Navier-Stokes  equations 

= 

transformed  components  of  Navier-Stokes  x-momentum  viscous 

flux 

Gb  Gn, 

G< 

= 

transformed  components  of  Navier-Stokes  y-momentum  viscous 

flux 

H<,  Hn, 

= 

transformed  components  of  Navier-Stokes  z-momentum  viscous 

flux 

hs 

= 

specie  enthalpy 
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1 

Iy 

J 

kB 

kfr*  kbr 

kP 

Keq 

li 

M 

N 

NC 

Nf 

Nr 

Ns 

P 

q 

R 

S 

Sa,  Sa-1 
Sb,  Sb1 
Sc,  Sc-1 

t 

T 

u 


=  identity  matrix 
=  intensity 

=  Jacobian  of  general  curvilinear  coordinate  transformation 
=  Boltzmann’s  constant  (Equation  28) 

=  forward  and  backward  reaction  rates,  respectively,  for  the  r*h 
reaction  path  (Equation  15) 

=  Planck's  constant  (Equation  2?) 

=  kfr/kbr  -  equilibrium  constant  for  the  r*h  reaction  (Equation  37) 
=  unit  vector  in  the  direction  of  Iv 
=  specie  molecular  weight 

=  total  number  of  harmonic  oscillators  (Equation  27) 

=  matrix  relating  coupled  Navier-Stokes/chemistry  conservative 
and  nonconservative  solution  vectors 
=  matrix  relating  fluid  dynamic  (Navier-Stokes)  conservative  and 
nonconservative  solution  vectors 
=  total  number  of  chemical  reaction  paths 
=  total  number  of  chemical  species 
=  pressure 
=  heat  flux 
=  universal  gas  constant 
=  entropy  (Equation  35) 

=  similarity  matrices  which  diagonalize  the  flux  Jacobian  A 

=  similarity  matrices  which  diagonalize  the  flux  Jacobian  B 

=  similarity  matrices  which  diagonalize  the  flux  Jacobian  C 

=  time 
=  temperature 

=  conservative  solution  vector  =  [p,  pu,  pv,  pw,  e]T 


')9 


U,  V,  w 

a,  o,  a 

V 

w 

x,  y,  z 

y+ 

z 


cartesian  components  of  velocity 
cartesian  components  of  specie  diffusion  velocities 
nonconservative  solution  vector  =  [p,  u,  v,  w,  ej]T 
chemistry  source  term  vector  =  [wi,  ...»  wNg,  0,  0,  0,  0]T 
cartesian  coordinate  system 
law-of-the-wall  coordinate 

curve-fitted  coefficients  for  the  thermodynamic  properties 
(Equations  33-36) 
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